Numerov matrix method with  SHO potential

(* Potential, desired max energy *)
V[s_]:=.5(s^2);ϵm=50.;
(* Determine grid *)
rturn=FindRoot[V[s]==ϵm,{s,ϵm}][[1,2]];d=(1/(Sqrt[2ϵm]));
n=Round[2((rturn/d)+4π )];s=Table[-((d (n+1))/(2 ))+d i,{i,n}];
(* Calculate KE matrix *)
𝕀[n_,d_]:=DiagonalMatrix[1+0Range[n-Abs[d]],d];
B=(1./12)(𝕀[n,-1]+10 𝕀[n,0]+𝕀[n,1]);
A=(1/( (d^2)))(𝕀[n,-1]-2 𝕀[n,0]+𝕀[n,1]);
KE=((-1)/2)Inverse[B].A;
(* Hamiltonian *)
H=KE+DiagonalMatrix[V[s]];
(* Energies, wavefunctions *)
{eval,evec}=Eigensystem[H];

(* Swap list ordering and show first 20 eigenvalues *)
in=Ordering[eval];

eval=eval[[in]];evec=evec[[in]];

eval[[;;20]]

(* Construct and plot analytical wavefunctions *)
psia[s_,n_]:=((π^(-1/4))/(Sqrt[(2^n)*(n)!]))Exp[-((s^2)/2)]*HermiteH[n,s];
analytical=Plot[-psia[s,9],{s,-10,10}];

Integrate[((psia[x,9])^2),{x,-∞,∞}]

(* Plot Numerov wavefunctions, normalize to 1 *)
f[i_]:=-(d (n+1))/2 +d i
scale=Map[f,Range[Length[evec[[10]]]]];
result=Transpose[{scale,evec[[10]]/(Sqrt[d])}];
numerov=ListPlot[result,PlotRange->All];

Show[numerov,analytical,ImageSize->500]