This Mathematica.SE question https://mathematica.stackexchange.com/q/284679/ is physical too, so I have decided to duplicate it here.
I have a set of anisotropic gaussian basis set which describes the ground state of the system with great accuracy.
The Hamiltonian of the system has the following form: $$H=-\frac{1}{2}\Delta-\frac{1}{r}+\frac{1}{2}\rho^2-1$$ where $r=(\rho,z,\phi)$ is a coordinate in the cylindrical system
The anisotropic gaussian basis set which I use to solve this task:$$\psi_j=e^{-b_{0j}\;z^2}e^{-a_{0j}\;\rho^2}$$where $a_{0j}$ and $b_{0j}$ are parameters
$a_{0j}$=0.500074, 0.502676, 0.510451, 0.507502, 0.52768, 0.625164, 0.935248, 1.826161, 3.845598, 8.293859, 18.014295, 39.20736, 85.3836, 185.975252, 405.095196, 882.399269, 1922.095421, 4186.829292, 9120.018288
$b_{0j}$=0.026284, 0.057253, 0.124713, 0.184064, 0.271658, 0.429143, 0.799226, 1.740925, 3.792194, 8.2604, 17.993331, 39.194226, 85.375371, 185.970096, 405.091965, 882.397245, 1922.094153, 4186.828497, 9120.01779
The exact value of the ground state energy from the literature is -1.0222124 and this basis set with great accuracy give this value (-1.02221). It's know the excited levels energies from the literature too: -0.174, -0.069, -0.036865, -0.022815, -0.015495. The basis set that I used give the follow values for the excited levels energies: -0.173374, -0.0199577, 0.312572, 1.18162... If the first excited level differs in the third decimal place, then the energies of the subsequent excited levels are completely different. This is because the basis set I am using is not complete and it is necessary to increase the number of basis functions in order to get more accurate values for the excited levels.
Is it possible, having a certain set of basic functions (in my case, it is 19 functions), to supplement it with additional functions so that it becomes more complete? How to do it? Because I would like to get the energies of the excited levels too.
The code in Wolfram Mathematica (in the code I rename $\rho\equiv r$):
ClearAll["Global`*"]
(the basis set)
b00 = {0.026284, 0.057253, 0.124713, 0.184064, 0.271658, 0.429143,
0.799226, 1.740925, 3.792194, 8.2604, 17.993331, 39.194226,
85.375371, 185.970096, 405.091965, 882.397245, 1922.094153,
4186.828497, 9120.01779};
a00 = {0.500074, 0.502676, 0.510451, 0.507502, 0.52768, 0.625164,
0.935248, 1.826161, 3.845598, 8.293859, 18.014295, 39.20736,
85.3836, 185.975252, 405.095196, 882.399269, 1922.095421,
4186.829292, 9120.018288};
b0 [j_] := b00[[j]];
a0 [j_] := a00[[j]];
Psi[r_, z_, j_] := Exp[-b0[j]z^2]Exp[-a0[j]*r^2];
nmax = Dimensions[b00][[1]];
Kk[r_, z_, j1_, j2_] :=
FullSimplify[
Psi[r, z, j2]*
Laplacian[Psi[r, z, j1], {r, [Theta], z}, "Cylindrical"]r2Pi];
Kx[j1_, j2_] := -(1/2)
NIntegrate[
Kk[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
Kxx = Table[Kx[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];
Ka[r_, z_, j1_, j2_] :=
FullSimplify[Psi[r, z, j2]1/2r^2Psi[r, z, j1]r2Pi];
KA[j1_, j2_] :=
KA[j1, j2] =
KA[j2, j1] =
NIntegrate[
Ka[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
KAx = Table[KA[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];
Ks[r_, z_, j1_, j2_] :=
FullSimplify[Psi[r, z, j2](-1)Psi[r, z, j1]r2*Pi];
KS[j1_, j2_] :=
KS[j1, j2] =
KS[j2, j1] =
NIntegrate[
Ks[r, z, j1, j2], {r, 0, Infinity}, {z, -Infinity, Infinity}];
KSx = Table[KS[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];
VP1[r_, z_] := -(1/Sqrt[r^2 + z^2]);
Px1[j1_, j2_] :=
Px1[j1, j2] =
Px1[j2, j1] =
NIntegrate[
Psi[r, z, j2]VP1[r, z]Psi[r, z, j1]r2*Pi, {r, 0,
Infinity}, {z, -Infinity, Infinity}];
Pxx1 = Table[Px1[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];
Bb[j1_, j2_] :=
Bb[j1, j2] =
Bb[j2, j1] =
NIntegrate[
Psi[r, z, j2]Psi[r, z, j1]r2Pi, {r, 0,
Infinity}, {z, -Infinity, Infinity}];
Bxx = Table[Bb[j1, j2], {j1, 1, nmax}, {j2, 1, nmax}];
EEE = Kxx + Pxx1 + KAx + KSx;
SortEigenvalues[{EEE, Bxx}]
Out[88]= {-1.02221, -0.173374, -0.0199577, 0.312572, 1.18162,
2.72596, 5.63718, 12.0333, 17.9735, 28.2709, 61.4758, 135.447,
298.338, 656.001, 1440.89, 3168.65, 7007.42, 15772.6, 37523.6}
