This Mathematica  code  is based  on  Trefethen' s  Program  8  from his book  Spectral  Methods  in  Matlab  but is extended using various refinements by GB that are used in the Phase Method

Warning: - spectral methods can be very effective for smooth potentials but are much less accurate for potentials that are discontinuous or even have low order discontinuous derivatives (e.g. breaks in slope). This is currently limited to Machine Precision (FP64). The Phase Method has no such limitations re smoothness, but it does require extended precision to get the same accuracy as spectral methods for smooth potentials.   I do not recommend simple finite difference methods for serious work. Finite Element Method with adaptive meshes can be quite effective however depending on implementation and choice of parameters.  The Phase Method basically just works.

$Assumptions={x ∈ (Reals)};

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]]])]

define potential and maximum energy level to calculate

k=2;ϵ=1;a=1;

U[x_]=ϵ/2 Abs[x/a]^k//FullSimplify;
etest=30;tpsf={2,2};(*turning point scale factors*)

The following cell estimates the number of bound states below etest, for homogeneous potentials of form 
U=ϵ/2 |x/a|^k with k>0.  
See Phase Method book, appendix A.  Ignore this for potentials that are not of this form.

η=ℏ/(a Sqrt[m ϵ]);μ=2k/(2+k);ℏ=1;m=1;
cc=.69;κ1=135;κ2=1427;(*from book appendix A*);
νest=1/2(cc Exp[-(2-k)/κ1]+(1-cc) Exp[-(2-k)/κ2]);
Ωest=Floor[((2 etest/ϵ)^(1/μ)(1/k)!)/((η Sqrt[π](1/μ)!))+νest];
Print[Framed["estimated number of bound states below ℰ="<>ToString[etest]<>" is "<>ToString[Ωest]]]

If[Ωest>200,Print[Ωest," is a whole lot of bound states. Quitting now."];Do[Beep[ ];Pause[.3],{5}];Quit[ ]];

This stuff attempts to determine appropriate range and sampling density for specified potential and upper limit etest.
It is modeled after things done in phase method:  setting an upper energy limit, finding the turning points at that energy,
choosing the range by scaling the turning points, and calculating the shortest wavelength in the range to determine sampling grid.
The phase method doesn’t actually require this, but it is used as an adequacy check for a simple version of  transition counting in the phase plot.

umin=FindMinimum[U[x],x]//Quiet//First

{xmin,xmax}=tpsf*MinMax[x/.NSolve[etest==U[x],x,Reals]]; λmin=2 π/Sqrt[2 (etest-umin)]; 
L=Ceiling[(xmax-xmin)/2]; δx=λmin/32.;   npts=2*Ceiling[L/δx];(*note:npts must be an even number*)
Print["L=",L," δx=",NumberForm[δx,3]," npts=",npts]

Limit[U'[x],x->0]

Plot[{etest,U[x]},{x,xmin,xmax},PlotRange->All]

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

k

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

The following tables are from Phase Method book, appendix F.

The spectral method agrees well with those values for k = positive even-integer e.g. 2,4,  
but much less precisely for other k values,  presumably because of the slope discontinuity at x=0 
from the Abs[x] function. The PM values are considered to be the accurate to about the number of digits shown.

for safekeeping these are iconized versions of the four tables below{,,,}

{{,},{,}}

Here are the wavefunctions of the 12 lowest energy states
Overall sign is immaterial.

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

(*comparison valid for harmonic oscillator only*)

Calculated eigenvalues and wavefunctions by any method will be inaccurate if they do not converge over the specified range. We mostly eliminate this concern by setting the range {xmin, xmax} based on an upper energy limit etest and the potential by scaling the  turning points as described in the Phase Method book.  The following plot shows the loss of accuracy of bound state energies at sufficiently high n value/

If[(U[x]===x^2/2)||U[x]===1/2Abs[x]^2,exactHO=Table[n-1/2,{n,1,Length[evals]}];
ListLogLogPlot[Abs[evals/exactHO-1],PlotLabel->"relative errors in bound state energy (harmonic oscillator potential)",
Frame->True,GridLines->Full,FrameLabel->{"n","rel error"}]]