define the function spectralsolve

RA[npts_]:=Range[1,npts];TR=Transpose;
spectralsolve[U_,npts_,L_]:=Module[{h,D2,col},(
h=2*Pi/npts; xray=h*RA[npts];xray=L*((xray-Pi)/Pi)//N;
col=Join[{-π^2/(3*h^2)-1/6},-(1./2)*(-1)^RA[npts-1]/Sin[h*RA[npts-1]/2]^2];
D2=(Pi/L)^2*ToeplitzMatrix[col];
{evals,evecs}=TR@SortBy[#[[1]]&]@TR[Eigensystem[-D2/2+DiagonalMatrix[U/@xray]]]
)]

test it out (results guaranteed to be pretty crummy for square well)

U[x_]:=If[Abs[x]<3,-10,0];
npts=2^14;L=16;
AbsoluteTiming[{evals,evecs}=spectralsolve[U,npts,L];]

MaxMemoryUsed[]/1024.^3 GB

scalefactor=Sqrt[Length[evals]/(2 L)];
wf[n_]:=Transpose[{xray,scalefactor*evecs[[n]]}]
GraphicsRow[{ListPlot[wf[1],PlotRange->{{-L,L},All},Joined->True],
ListPlot[wf[9],PlotRange->{{-L,L},All},Joined->True]}]

npts

select energies states below threshold ℰmax

ℰmax=0;

bs=Select[evals,#<ℰmax&]

Iconize[bs]

now compare with Phase Method eigs WP25 and essentially exact rootfinding (calculated WP 100)

PMeigsWP16=(*PM WP 16*)

PMeigsWP25=(*PM WP 30*)

RFeigs=(*rootfinding WP100*)

Abs[bs/RFeigs-1]

Abs[PMeigsWP16/RFeigs-1]

ListLogPlot[{Abs[bs/RFeigs-1],Abs[PMeigsWP16/RFeigs-1],Abs[PMeigsWP25/RFeigs-1]},PlotRange->All,Frame->True,GridLines->Full,PlotMarkers->Automatic,
PlotLegends->{"relative error spectral","relative error PM WP16","relative error PM WP25"},ImageSize->10*72]