Why occur a giant increasing of absolute value of minimal eigenvalue with rise of the basis functions number?

Question

I have the following code:

code 1

In[193]:= ClearAll["Global`*"]

nmax = 5;

(*3d oscillator basis (l=0)*)
Psi[r_, n_] := (-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n + 3/2]]* 
   LaguerreL[n, 1/2, r^2];
(*kinetic energy*)
Kk[r_, n1_, n2_] := 
  FullSimplify[
   Psi[r, n2]*
    Laplacian[Psi[r, n1], {r, \[Theta], \[Phi]}, "Spherical"]];
Kx[n1_, n2_] := -1/2*Integrate[Kk[r, n1, n2]*r^2, {r, 0, \[Infinity]}];
KK = Table[Kx[n1, n2], {n1, 0, nmax}, {n2, 0, nmax}];

(*potential energy*)
VH1[r_] := -1/r;
VH2[r_] := -(1/(2*r))*Exp[-r*1.5];
Px1[n1_, n2_] := 
  Integrate[Psi[r, n2]*VH1[r]*Psi[r, n1]*r^2, {r, 0, \[Infinity]}];
Px2[n1_, n2_] := 
  Integrate[
   Psi[r, n2]*FullSimplify[VH2[r]]*Psi[r, n1]*r^2, {r, 
    0, \[Infinity]}];
Px[n1_, n2_] := Px1[n1, n2] + Px2[n1, n2];

PP = Table[Px[n1, n2], {n1, 0, nmax}, {n2, 0, nmax}];
EE = Min[Eigenvalues[KK + PP]]

Out[205]= -0.675029

But Mathematica calculates the code 1 slowly, therefore I have rewritten it and receive the code 2. Mathematica calculates the code 2 much faster so I can use a lot of basis functions.

If I use 25 basis functions (nmax=24) EE = -0.712228, but if I put nmax=25 value of EE become -13.236, and this value is unrealistic because from general theory we know that with increase of basis functions number the minimal eigenvalue should be change less and less.

Why occur such giant increasing of absolute value of EE? What wrong in the code2?

code 2

In[207]:= ClearAll["Global`*"]

nmax = 5;

(*3d oscillator basis*)
(*Psi[r_,n_]:=(-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n+3/2]]* \
LaguerreL[n,1/2,r^2];*)

(*kinetic energy*)
nmx = nmax;
hm = ConstantArray[0, {nmx + 1, nmx + 1}];
Do[hm[[1 + n1, 1 + n2]] = 
  KroneckerDelta[n1 - n2]*(n1 + 3/4) - 
   Sqrt[(2 (n1)^2 + (n1))/8] KroneckerDelta[n1 - n2 - 1];
 hm[[1 + n2, 1 + n1]] = hm[[1 + n1, 1 + n2]];, {n1, 0, nmx}, {n2, 0, 
  n1}]
KK = N[hm];

(*potential energy*)

f0[a_, n_] = 
  Integrate[r^n*Exp[-r^2]*Exp[-a*r], {r, 0, Infinity}, 
   Assumptions -> {a > 0, n > 0}];
f[a_, n_] := f[a, n] = N[f0[a, n], 100];
coeff[n_] := (-1)^n*Sqrt[2*n!/Gamma[n + 3/2]];
(*Coulomb part VH1[r]*)
Px1[n1_, n2_] := 
  Px1[n1, n2] = 
   If[n1 > n2, Px1[n2, n1], 
    coeff[n1]*coeff[n2]*(-1)*
      Total[Map[Function[{x}, f[0, x[[1, 1]] + 1]*x[[2]]], 
        CoefficientRules[
         LaguerreL[n1, 1/2, r^2]*LaguerreL[n2, 1/2, r^2], r]]] // N];

(*VH2[r]*)
Px2[n1_, n2_] := 
  Px2[n1, n2] = 
   If[n1 > n2, Px2[n2, n1], 
    coeff[n1]*coeff[n2]*(-(1/2))*
      Total[Map[Function[{x}, f[-1.5, x[[1, 1]] + 1]*x[[2]]], 
        CoefficientRules[
         LaguerreL[n1, 1/2, r^2]*LaguerreL[n2, 1/2, r^2], r]]] // N];

Px[n1_, n2_] := Px1[n1, n2] + Px2[n1, n2];
PP = Table[Re[Px[n1, n2]], {n1, 0, nmax}, {n2, 0, nmax}];
EE = Min[Eigenvalues[KK + PP]]

Out[220]= -0.675029