This is Mathematica  code  based  on  Trefethen' s  Program  8  
from Spectral  Methods  in  Matlab

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 the code for harmonic oscillator

note : npts  must  be  an  even  number

U[x_]:=x^2/2;npts=2^8;L=16;
{evals,evecs}=spectralsolve[U,npts,L];
scalefactor=Sqrt[Length[evals]/(2 L)];
wf[n_]:=Transpose[{xray,scalefactor*evecs[[n]]}]

TableForm@Partition[Take[NumberForm[#,{16,12}]&/@evals,150],5]

here are the lowest 12 wavefunctions -

GraphicsGrid[Partition[Table[ListPlot[wf[n],PlotRange->{{-L,L},All},Joined->True],{n,1,12}],4],ImageSize->10*72]

in discretization (e.g. spectral) methods
for large enough n eigenvalues and wavefunctions quickly lose accuracy
if the nth wavefunction isn’t fully converged to near zero over the x-range
this criterion can be used as a guide to assessing the valid range of n for eigenvalues and wavefunctions

exactHO=Table[n-1/2,{n,1,npts}];

ListLogLogPlot[Abs[evals/exactHO-1],PlotLabel->"relative errors (harmonic oscillator potential)",Frame->True,GridLines->Full,FrameLabel->{"n","rel error"}](*comparison valid for harmonic oscillator only*)