I've been working on a toy model of a variational method of a Hamiltonian H with an overlap matrix S. The point is that when H and S are small enough, Mathematica gives the correct eigenvalues (≈ -1,-0.25,-0.111... hydrogen-like eigenvalues).
The problem is that when i increase the basis, the correct eigenvalues become zero. For example, if my basis is made of 34 states i have some correct eigenvalues of the exact problem and some others (positive ones) that are wrong, but are expected, however, if i increase the basis to 35 states a few of the correct eigenvalues become 0, and if the basis is, let's say 50 states, all the correct eigenvalues become 0.
I know i don't need such amount of states for this problem but the real problem i want to address is more complex than this one and i'd like to understand what is happening here.
Here i put the code i am using, S[i,j] and H[i,j] are the matrix elements of Sij[n] and Hij[n], where n is the size of the basis:
S[i_, j_] :=8 Sqrt[1.7^(3 i + 3 j - 6)]/(1.7^(i - 1) + 1.7^(j - 1))^3
H[i_, j_] := -8 Sqrt[1.7^(3 i + 3 j - 6)]/(1.7^(i - 1) + 1.7^(j -
1))^3*(0.01*(1.7^(i - 1) + 1.7^(j - 1)) - 0.01^2*1.7^(i + j - 2))
Sij[n_] := Table[S[i, j], {i, 1, n}, {j, 1, n}]
Hij[n_] := Table[H[i, j], {i, 1, n}, {j, 1, n}]
Eig[n_]:=Sort[Eigenvalues[{Hij[n],Sij[n]}]
I have tried not to use directly the Mathematica default Eigenvalues algorithm, for example using the inverted Sij[n] or the Cholesky decomposition, but i had similar results.