(*-------------------------------
this code is a little more complex than the python example 
since it does arbitrary precision, parallel, more graphics, and tests 
and much of it is comments like these. 
This example code does not do wave functions however, just energies. 
If pasted into Mathematica,  best put code chunks into separate cells. 
This code is formatted for inclusion in portrait mode document but 
Mathematica Notebook format is much preferred.  
-------------------------------*)

(*initial setup; choose desired digits of working precision <42 or so;  
$MachinePrecision is one option. More digits slows the computation 
but is usually more accurate. Too high a setting however can make 
eigenvalue search twitchy (unstable)*)
workingprecision=30;

(*estimated ultimate accuracy, based on experience; 
your mileage may vary*)
ngooddigits=workingprecision/2; 

(*set nominal value and slope initial conditions at xmin
these are exact numbers so precision isn't degraded
otherwise use SetPrecision*)
Phi0=1/10;Phi0p=1/10; 

(*define the potential here however you want to. 
Here we just define a finite square well using an If statement*)
a=5;V0=-20;
U=If[Abs[x/a]<1,V0,0]; 

(*set etest to the upper end of your energy search range*)
etest=-1/10;

(*set the range - should be well outside turning points at etest. 
nearly-unbound states will require wider x-range for accuracy*)
{xmin,xmax}={-40,40}; 

(*verify thatg,xmax are OK by checking the phaseplot at etest*)
{umin,uminxsoln}=NMinimize[U,x,Method->"RandomSearch"]//Chop;
uminx=x/.uminxsoln;
Print["here is the potential plot:"]
Plot[{U, etest},{x,xmin,xmax},PlotRange->All,Frame->True]

(*set energy search range for eigenvalue search. 
set to whatever makes sense for your problem*)
{eminsearch,emaxsearch}={umin,etest};
searchinterval=SetPrecision[{eminsearch,emaxsearch},workingprecision];

(*the "nbs"  (aka "omega") function estimates the number of bound states below 
specified energy; it also calculates the phaseplot in the process. 
clear out previously memo-ized values*)
ClearAll[nbs];

(*sampling grid for counting transitions, just has to be small enough. 
This could be estimated from max momentum at potential minimum, but no need.*)
dx=.001;

(*make the variables extended-precision so as to not degrade precision.
A nicer way to count transitions is to use "WhenEvents" in NDSolve, 
but the approach shown here is portable to other languages.*)
{Phi0,Phi0p,xmin,xmax,e}=SetPrecision[{Phi0,Phi0p,xmin,xmax,e},workingprecision];
nbs[e_]:=nbs[e]=(soln=Quiet@NDSolve[{-1/2Phi''[x]+(U-e)Phi[x]==0,
Phi'[xmin]==Phi0p,Phi[xmin]==Phi0},Phi[x],{x,xmin,xmax},MaxSteps->10^7,
Method->"LSODA",WorkingPrecision->workingprecision][[1]];
phaseplot[x_]=Tanh[Phi[x]/.soln];
transitions=DeleteCases[0]@Differences@Quiet@Table[Sign[Tanh[Phi[x]/.soln]],
{x,xmin,xmax,dx}]/2;
Return[{e,nboundstates=Length[transitions]}]);

(*examine the phaseplot before running the rest of calculation 
to verify the x-range is suitable for the chosen maximum energy. 
transitions should be sharp, not rounded*)
Print["here is the Phase Plot at upper energy limit etest=",etest]
nbs[etest]; 
dophaseplot:=Plot[phaseplot[x],{x,xmin,xmax}, Frame->True,
PlotLabel->"# bound states="<>ToString[nboundstates],
PlotRange->All,PlotPoints->2000,ImageSize->8*72];

(*let's do a test run to make sure the x-range is set sensibly. 
if transitions are rounded instead of sharp, extend xmin and xmax 
past the turning points; if overflow is encountered, reduce the range 
or better, choose xmin by scaling the turning point values
or scale down the initial slope and value boundary conditions*)

dophaseplot

(*define "sift" and "refine" functions
accuracy seems to be limited to about half the working precision*)
maxiterations=Min[100,Ceiling[(workingprecision/2*Log[2,10]+5)]+2];
(*"sift" maps out where eigenvalues are located (serial process)*)
sift=bsearch/@(Partition[Flatten[#],2])&;
(*"refine" gets more accurate estimates (parallel process)*)
refine=FixedPoint[sift,#,maxiterations]&;

(*define evolutionary tree-like binary search function*)
ns[e_]:=nbs[e][[2]];intervals={};
bsearch[intrvl_]:=Module[{vlow,vhigh,vmid,mid,newintervals},(
{elow,ehigh}=intrvl;emid=(elow+ehigh)/2;
{vlow,vmid,vhigh}={ns[elow],ns[emid],ns[ehigh]};newintervals={};
If[vlow<vmid,AppendTo[newintervals,{elow,emid}]];
If[vmid<vhigh,AppendTo[newintervals,{emid,ehigh}]];
Return[newintervals]);]

(*first, sift for eigenvalues (serially), enough to resolve the eigenvalues. 
this is mostly for visualization purposes.  Don't worry if sifting is incomplete -- 
the refine process (in parallel) will still take care of the rest*)

Print["sifting can take a few minutes, please be patient..."];
(*rough estimate of number of bisections needed to separate 
the eigenvalue intervals*)
snit=Ceiling[Log[2,nboundstates]]+5;
{sifttime,intervals}=
AbsoluteTiming[Flatten[Nest[sift,searchinterval,snit],1]];
Print["sifting complete: sifting time=",sifttime," sec,  # of intervals=",
Length[intervals]];
ListPlot[
Interval/@intervals,PlotRange->All,Frame->True]

(*this optional plot shows the interval search, using values already "memoized"*)
Print["the following plot shows the evolution of the binary search process"]
siftplot=ListPlot[Interval/@Partition[Flatten[NestList[sift,searchinterval,10]],2],
PlotRange->All,ImageSize->Large,Frame->True,
FrameLabel->{{Subscript[\[ScriptCapitalE],n],None},
{None,"U="<>ToString@TraditionalForm[U]}}]

(*then "refine" the found eigenvalues*)
Print["refining also can take a few minutes, please continue to be patient..."];
{refinetime,intervals}=AbsoluteTiming[
Flatten[ParallelMap[refine,intervals]//Quiet,2]];
eigs=Mean/@intervals;neigs=Length[eigs];roundedeigs=N[eigs,ngooddigits];
Print["sift time=",sifttime," refine time=", refinetime," total=",
sifttime+refinetime," (seconds)"];
Print["neigs=",neigs,"  ",(sifttime+refinetime)/neigs," seconds/eigenvalue"];

(*here are the eigenvalue estimates which are 
typically accurate to about workingprecision/2 digits
your mileage may vary*)
TableForm[Transpose[{Range[neigs],roundedeigs}]]

(*--------------------
the following code checks accuracy of the calculated eigenvalues 
for a symmetric finite square well with the potential parameters above. 
the following code can be ignored or deleted for other potentials
--------------------*)
setpars={m->1,hbar->1};
setvars={k->Sqrt[2 m (e-V0)]/hbar,Gamma->Sqrt[2 m (-e)]/hbar};
zeta=hbar^2/(2 m a^2);
nbound[V0_]:=1+Floor[(2/Pi)Sqrt[-V0/zeta]]/.setpars;
Print["nbound states=",nbound[V0]];

(*applying the following function to an eigenvalue should give
zero for even parity states*)
eventest[e_]=Cot[k a]-k/Gamma/.setvars/.setpars;
(*applying the following function to an eigenvalue should give 
zero for odd parity states*)
oddtest[e_]=Tan[k a]+k/Gamma/.setvars/.setpars;
evens=eventest[eigs[[1;;;;2]]];
odds=oddtest[eigs[[2;;;;2]]];
Print["est number of good digits=",ngooddigits];
Print["this is a plot of the absolute error in the eigenvalues"];
ListLogPlot[Tooltip[Abs@Riffle[evens,odds]], Frame->True]
(*these values should be close to zero if the calculated 
eigenvalues are correct*)
(*-------------------------------------------*)