Hydrogen Atom using Numerov Method

Due to (⟨r⟩~(n^2), we need a variable mesh grid.  We then need a corresponding rescaling of the wavefunction to preserve the Numerov form.  Let ψ(r)=J(s(r))χ(s(r)).  We want to choose J so that first derivatives are eliminated in the Schrodinger equation.

One choice of mesh is r=(s^2)

(∂_(r,r))(J[(Sqrt[r])]χ[(Sqrt[r])])/.{r->(s^2)}//PowerExpand

Collect[%,χ'[s]]

To eliminate the first derivative we need 2 J[s]+8 s J'[s]=0

DSolve[-((J[s])/(4 (s^3)))+((J'[s])/(2 (s^2)))==0,J[s],s]

To get in eigenvalue form, multiply SE by 1/(Sqrt[s])

SE=(1/(Sqrt[s]))(((-1)/2)(∂_(r,r))((r^(1/4))χ[(Sqrt[r])])-(ℰ-V[r])(r^(1/4))χ[(Sqrt[r])]/.{r->(s^2)})//PowerExpand//Simplify//Expand

V[r_]:=((-1)/r);ϵm=20.;
(*determine grid, rmax=4000*)
n=500;ds=(Sqrt[4000.])/n;
s=Table[-ds /2+ds i,{i,n}];
(*calculate KE matrix*)
one[n_,d_]:=
DiagonalMatrix[1+0Range[n-Abs[d]],d];
B=(one[n,-1]+10 one[n,0]+one[n,1])/12.;
A=(one[n,-1]-2 one[n,0]+one[n,1])/(ds^2);
KE=-DiagonalMatrix[(1/(8(s^2)))].Inverse[B].A;
(*Hamiltonian*)
H=KE+DiagonalMatrix[V[(s^2)]+(3/(32 (s^4)))];
(*Energies, wavefunctions*)
{eval,evec}=Eigensystem[H];in=Ordering[eval];eval=eval[[in]];evec=evec[[in]];

calculate effective quantum number from the energies, accurate out to about n=40

(Sqrt[((-1)/(2eval[[;;50]]))])

Here is the n=40 wavefunction as a function of r

ListPlot[{(s^2),(s^(1/2))evec[[40]]}\[Transpose],PlotRange->All,Joined->True]

Plotted vs s, the nonlinear grid has smoothed out the wavefunction to be closer to a sine wave

ListPlot[{s,evec[[40]]}\[Transpose],PlotRange->All,Joined->True]

1s wavefunction still pretty accurate

Show[ListPlot[{(s^2),(s^(1/2))evec[[1]]}\[Transpose],PlotRange->{{0,6},All},Joined->False],Plot[x Exp[-x],{x,0,6},PlotStyle->Red]]