Why do different calculations give different eigenvalues?

Question

I find the eigenvalues of the Hamiltonian in two ways - numerically and matrix . Why do different calculations give different eigenvalues?
The Hamiltonian is H = -1/2*Laplacian -(2/r)
It can be seen that the minimum energy calculated numerically is -2.00001, and matrix-wise is -1.55846. Why is this happening?

1. numerically

In[82]:= H = 
  Simplify[-1/2*
     Laplacian[Psi[r], {r, \[Theta], \[Phi]}, "Spherical"] - 
    Psi[r]*2/r];

{vals, funs} = 
  NDEigensystem[{H + Psi[r]*0.5}, Psi[r], {r, 0, 100}, 100, 
   Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.015}}}, 
     "Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}];

In[84]:= (Take[Sort[vals], 5] - 0.5)

Out[84]= {-2.00001, -0.500001, -0.222223, -0.125, -0.0800001}

2. matrix

In[26]:= ClearAll["Global`*"]
Psi1[r_, n1_] = (2 E^(-(r/n1)) Sqrt[n1!/(-1 + n1)!] Hypergeometric1F1[
      1 - n1, 2, (2 r)/n1])/n1^2;

Psi2[r_, n2_] = (2 E^(-(r/n2)) Sqrt[n2!/(-1 + n2)!] Hypergeometric1F1[
      1 - n2, 2, (2 r)/n2])/n2^2;

(*kinetic energy*)

K[r_, n1_, n2_] = 
  FullSimplify[
   Psi2[r, n2]*
    Laplacian[Psi1[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];

(*potential energy*)

P[r_] = -2/r;

(*calculation of matrix elements,I am using 20 basis functions*)

EE = Table[-1/2*
      NIntegrate[K[r, n1, n2]*r^2, {r, 0, \[Infinity]}, 
       MaxRecursion -> 15, WorkingPrecision -> 32] + 
     NIntegrate[Psi2[r, n2]*P[r]*Psi1[r, n1]*r^2, {r, 0, \[Infinity]},
       MaxRecursion -> 15, WorkingPrecision -> 32], {n1, 1, 20}, {n2, 
     1, 20}] // Chop;

(*eigenvalues*)

Eeig = Eigenvalues[EE]

Out[32]= {-1.5584632981876349185146804146007, \
-0.3623635914395303877441420153245, \
-0.15669498305937198905543539839599, \
-0.08683775701950283493030149036754, \
-0.05501734483287334369748266587302, \
-0.03790244105929780721020220539079, \
-0.027651029147976237708049044437355, \
-0.021028273261568554013526089309235, \
-0.016502425356337140815385665090984, \
-0.013271998387649143824131249079151, \
-0.010884520722145291443991307150682, \
-0.009068777306934890312009174759412, \
-0.007654074742951805746421678011622, \
-0.006528495175126375499971622294854, \
-0.005615980338382413802381791780753, \
-0.004863061581404067179392515633070, \
-0.004230715950120777335031723062825, \
-0.003688843483376281574694297474038, \
-0.0032113000791884747824869858880628, \
-0.0027659547381622397716009741638790}

In[33]:= 
Emin = Min[Eeig]

Out[33]= -1.5584632981876349185146804146007

Answer 1

A common trick in spherical coordinates is to define $\psi(r)=\frac{u(r)}{r}$ to simplify the normalization and Laplacian expressed in terms of $u(r)$:

$$ \nabla^2\psi(r)=\psi''(r)+\frac{2}{r}\psi'(r) \qquad\leftrightarrow\qquad \nabla^2\frac{u(r)}{r}=\frac{u''(r)}{r}\\ \int_0^{\infty}\psi^2(r)r^2dr=1 \qquad\leftrightarrow\qquad \int_0^{\infty}u^2(r)dr=1 $$

The Schrödinger equation then becomes $$ -\frac12\nabla^2\psi(r)-\frac{2}{r}\psi(r)=E\psi(r) \qquad\leftrightarrow\qquad -\frac12u''(r)-\frac{2}{r}u(r)=E u(r) $$ which is much easier to deal with (looking like pseudo-straight coordinates with trivial Jacobian).

H = -1/2 u''[r] - 2/r u[r];

{vals, funs} = NDEigensystem[H, u[r], {r, 0, 100}, 100,
  Method ->
    {"SpatialDiscretization" -> {"FiniteElement",
        {"MeshOptions" -> {"MaxCellMeasure" -> 0.015}}},
     "Eigensystem" -> {"Arnoldi", "MaxIterations" -> 10000}}] //
  Transpose // Sort // Transpose;

vals[[;; 5]]
(*    {-1.46467, -0.424888, -0.199048, -0.11502, -0.0748264}    *)