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}    *)