Since A is non-singular, we can write (warning: below is math, not code)
M = A + s*D
=> M/A = I + sD/A (where D/A stands for D*inv(A))
=> lambda*M/A = lambda*I - D/A (where lambda = -1/s)
Therefore you need to solve the eigenvalue problem:
for lambda. You can find lambda using:
setting s = -1./lambda, you have n solutions (with multiplicity) for d = 0, where n is the size of A.
NOTE: Don't try to find eigen vectors (you dont need them here anyway) using Matlab's eig(), since it wont work if eigenvalues are repeated.
Now for non-zero d, we want:
d/det(A) = det(I + s*D/A)
since det(D/A) = product(lambda), we have to choose s so that
product_over_i(1 + s*lambda(i)) = d/det(A)
So your "short-cut" procedure is:
1. Find vector of eigenvalues, lambda, of D/A
2. Solve prod(1 + s*lambda) = d/det(A)
(1) should not be a problem even for very large matrices.
I'm not sure about (2) unless n is not too large. You could try the roots() or fzero() functions in Matlab. Note that if you use roots(), you will have to compute the coefficients of the polynomial first from the product in (2) and you may lose a lot of precision in computing the coefficients themselves. I think it would be better to use fzero().
Also note that there are n solutions in general, so you'll have to find a way to choose the solution based on your needs and also make an appropriate initial guess.
[edited to correct sign error].