For comparison this notebook uses Mathematica’s NDEigensystem to solve the problem, using finite element method.

This notebook also uses some strategies employed in the phase method such as setting an upper energy limit and 
setting the range by scaling classical turning points. It also by necessity sorts the levels because they arrive out of order,
and calculates an excess number of levels so as to get the relevant ones.  It seems to work OK if the number of states needed isn’t large.

(*U[x_]:=x^2/2;*)
U[x_]:=If[Abs[x]<3,-10,0];
H=-1/2 D[#,{x,2}]+U[x]* # &;
etest=0;(*upper limit for energy levels sought*)

Clear[xmin,xmax];tpsf={2,2};
Quiet[{xmin,xmax}=tpsf*(x/.NSolve[U[x]==etest,x])];
If[!And@@(NumericQ/@{xmin,xmax}),Print["couldn't calculate turning points"];{xmin,xmax}={-10,10}];
Print["xmin,xmax=",{xmin,xmax}]

The following code generates a bunch (ntries) at various energiess, and collects states those with eigenvalues below etest.
For negative energy states, to catch all of the relevant ones,  you may have to calculate several times the number you are actually seeking, and filter out the irrelevant ones.

ntries=40;(*number of states to calculate - sometimes this must exceed the number of states you actuall seek*)
rawevals=.;
{evals,evecs}=NDEigensystem[H[ψ[x]],ψ[x],{x,xmin,xmax},ntries,Method->{"PDEDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->.001}}}];
rawevals=evals;
thestates=Select[SortBy[Transpose[{evals,evecs}],#[[1]]&],#[[1]]<etest&];
nstates=Length[thestates];
Print["nstates found below etest=",nstates];
{evals,evecs}=Transpose[thestates];
Print["evals=",evals];
Print[NumberLinePlot[rawevals]];
Plot[evecs[[1;;3]],{x,-10,10},PlotRange->All]

orthogonality table

ParallelTable[NIntegrate[evecs[[i]]evecs[[j]],{x,xmin,xmax}]//Quiet,{i,1,nstates},{j,1,i}]//Chop[#,10^-7]&//TableForm