The problem is expressed as A*X = B*X*Lambda, where both A and B are
real symmetric (or complex Hermitian) and B is positive-definite. Lambda
is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
After solve, A is overwritten with Lambda, and B is overwritten with I.
No intermediate state is implemented, nor compact storage.