Phase Method basic eigenvalue search

Grant B. Bunker, Illinois Institute of Technology Feb 20, 2026
See book “The Phase Method” Kindle Direct Publishing 2004
To try this, just choose “Evaluation/Evaluate Notebook” to run it.  Typical run time is on the order of minutes for 10 or 20 eigenvalues, accurate to 30 decimal places. 
You may want to do Evaluation/Quit Kernel and Cell/Delete All Output to start with a clean slate first. 
Modified from example code in book “The Phase Method”.
ODE solver changed from LSODA to StiffnessSwitching.
Transition finding is done within ODE solver itself using WhenEvents.
If no minimum to the potential function (e.g. Coulomb, Yukawa) is found, a rough estimate is generated to give a lower limit to the energy search range.

This code, while still very basic, is a little more complex than the python example in the book since it does:
arbitrary precision, parallel computation, better graphics, and some accuracy tests, and much of it is comments like these. 
However this example code does not do wave functions, just energies. 
However a “workhorse” notebook does do that, and a lot more, maybe too much — as of Nov 29, 2025, this can be downloaded at https://gbxafs.iit.edu/phase-method/
The notebook that you are now reading also does not attempt to automatically estimate a good x-range by scaling the extremal classical turning points,  as does the “workhorse notebook”.  
Setting the x-range appropriately is up to you, the user, to do explicitly — but the phase plot can be instructive for doing that. 
A progress pane is shown during the parallel refinement to show you how it’s going.  

There is code at the end that is only appropriate for comparison with precise calculations for a finite square well as described below. 
[Feb 20 2026 - small improvements made to the code comparing PM results with “exactish” rootfinding results for finite square well at end of notebook.]

initial setup

ClearAll["Global`*"]

starttime=AbsoluteTime[ ];

speakelapsedtime=True;

Choose desired digits of working precision $MachinePrecision ≤ workingprecision ≤ 150 or so;
$MachinePrecision is one option (fastest). More digits make the computation a bit slower but far more accurate.
Normally, higher accuracy requires extending the x-range, especially for levels near the continuum.

workingprecision=60;(*make this twice the desired number of accurate digits; a multiple of 3 makes for nicer format tables*)
integrationmethod="StiffnessSwitching";

estimated ultimate accuracy, based on experience;  your mileage may vary

ngooddigits=Floor[workingprecision/2];

set nominal value and slope initial conditions at xmin
for singular potentials e.g. Coulomb choose Φ0=0, Φ0p=10^4 or so.

Φ0=1/10;Φ0p=1/10;

Define the potential here however you want to.  Here as an example we just define a finite 1D square well using an If statement.
Add l(l+1)/(2 x^2) for higher l values in (say) 3D but the test code at the end isn’t appropriate for that case. 
In Mathematica you can use rational fractions to retain (infinite) precision, or use SetPrecision as used below.

γ=1;(*inverse mass factor -- normally set this to 1*)

a=3;V0=-10;(*half-width and depth*)
U=If[Abs[x/a]<1,V0,0];

Ustring=ToString[TraditionalForm@U]

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.

WP=workingprecision;

{xmin,xmax}={-30,30};
{xmin,xmax}//N

you should verify that xmin, xmax are OK by checking the phaseplot at etest
this code attempts to find minimum in potential to estimate search range;  
if there is no minimum, e.g. -1/x potential,  it uses feynman’s heuristic to estimate one. 
For such potentials set xmin to be a small but finite number e.g. 10^(-workingprecision/2) and Φ0=0 there, and Φ0p=10^4 or so.

{umin,uminxsoln}=Check[NMinimize[U,x,Method->"RandomSearch"]//Chop,NMinimize[1/(2x^2)+U,x],NMinimize::cvdiv];
Print["here is the potential plot and etest - the turning points are where they intersect:"]
Plot[{U,etest},{x,xmin,xmax},PlotRange->All,Frame->True]

Set energy search range for eigenvalues to whatever makes sense for your problem.
Set precision of these variables to workingprecision so as not to degrade precision of computation.

{eminsearch,emaxsearch}={umin,etest};
searchinterval=SetPrecision[{eminsearch,emaxsearch},workingprecision];

the "nbs" (number of bound states aka "Ω") function calculates the number of bound states below specified energy; it also generates the phase plot in the process. 
Clear out previously cached (“memo-ized”) values if they exist.

ClearAll[nbs];

Make the variables extended-precision so as to not degrade precision in the calculation
note: all of the other input parameters need similar care!

{Φ0,Φ0p,xmin,xmax,U}

{Φ0,Φ0p,xmin,xmax,U}=SetPrecision[{Φ0,Φ0p,xmin,xmax,U},workingprecision];

nbs[e_]:=nbs[e]=(Ω=0;soln=Quiet@NDSolve[{-γ/2 Φ''[x]+ (U-e)Φ[x]==0,Φ[xmin]==Φ0,Φ'[xmin]==Φ0p,WhenEvent[Φ[x]==0,Ω++]},Φ[x],{x,xmin,xmax},MaxSteps->10^7,Method->integrationmethod,WorkingPrecision->workingprecision][[1]];
phaseplot[x_]=Tanh[Φ[x]/.soln];
Return[{e,Ω}]);

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 and there should extended regions of “saturation” at the intervals by xmin and xmax.

nbs[etest];
dophaseplot:=Plot[phaseplot[x],{x,xmin,xmax},Frame->True,PlotLabel->"# bound states="<>ToString[Ω],PlotRange->All,PlotPoints->2000,ImageSize->8*72];

sanity check

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 further past the turning points; 
If overflow is encountered, reduce the range, or better, choose xmin and xmax by scaling the turning point values as described in the book, “The Phase Method”.

Print[Framed["Here is the phase plot at the upper energy limit etest="<>ToString[TraditionalForm@etest]]];
dophaseplot

start search

define "sift" and "refine" functions — accuracy is limited to about half the working precision so choose that appropriately

maxiterations=Max[100,Ceiling[(workingprecision/2*Log[2,10]+5)]]

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

sift (serial process, as implemented)

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[Framed@"for some potentials, 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,Ω]]+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[ℰ,n],None},{None,"U="<>Ustring}}]

refine (parallel process)

Then "refine" the found eigenvalues

Print[Framed@"for some potentials, refining also can take a few minutes\n ... maybe enough time to get a cup of coffee or tea"];
{refinetime,intervals}=AbsoluteTiming[Flatten[ParallelMap[refine,intervals]//Quiet,2]];
eigs=Mean/@intervals;neigs=Length[eigs];roundedeigs=N[eigs,ngooddigits+2];
Print["sifting time=",sifttime,"; refining time=",refinetime,"; total=",sifttime+refinetime," (seconds)"];
Print["neigs=",neigs,";  ",(sifttime+refinetime)/neigs," seconds/eigenvalue"];

Here are the eigenvalue estimates, which are pretty reliably accurate to about workingprecision/2 digits as long as the x-range is sufficient

neigs

ngd=ngooddigits;NumberForm[TableForm[Transpose[{Range[1,neigs],roundedeigs}]],{ngd,ngd},DigitBlock->3]

The following iconized object contains eigenvalues rounded to correct precision. 
It can be copied into other Mathematica notebooks for further computation.

Iconize[roundedeigs]

elapsedtime=Round[(AbsoluteTime[]-starttime)/60.,.1] min;

Print[Framed@("total elapsed time= "<>ToString[elapsedtime])];
{$Version,$OperatingSystem,$ProcessorType,$ProcessorCount};

Now[] (*time stamp*)

delete this audio output if you don’t want it, or set the flag speakelapsedtime=False at the start of this notebook

speakthis="total elapsed time was "<>ToString[elapsedtime];
If[speakelapsedtime,AudioPlay[SpeechSynthesize[speakthis]];]

Optional: the following code can be ignored or deleted for other potentials. 
As a verification test, the following code checks accuracy of the calculated eigenvalues 
for a symmetric 1D finite square well with the potential parameters above,
equations standard QM methods,  from GB PHYS 509 notes. 
The following does not account for the inverse mass parameter γ if that is used. This Gamma is not that γ.

opts={Frame->True,FrameLabel->{"n","ℰ"},GridLines->Full,ImageSize->6*72};

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;

solve the exact equations for eigenvalue by numerical rootfinding, 
here using PM determined eigenvalues with random perturbations as starting guesses
and scale determined by lowest eigenvalue splitting

rand:=(eigs[[2]]-eigs[[1]])*.1*Random[]

eroots=e/.FindRoot[eventest[e],{e,#},WorkingPrecision->workingprecision]&/@(eigs[[1;;;;2]]+rand);
oroots=e/.FindRoot[oddtest[e],{e,#},WorkingPrecision->workingprecision]&/@(eigs[[2;;;;2]]+rand);
exactish=Riffle[eroots,oroots];

Precision/@exactish(*precision of "exact" benchmark values*)

these are the maximum and minimum relative errors of the PM-computed eigenvalues relative to the “exactish” rootfinding solutions to the exact equations

ngooddigits

{minerr,maxerr}=MinMax[Abs[eigs-exactish]];

If[maxerr<10^(-ngooddigits),Print[Framed@Framed@"ngooddigits criterion is OK"],
Print[Framed@Framed@"ngooddigits criterion not met - please check your x-range"]]

ngooddigits

N[{minerr,maxerr},2](*best and worst accuracy of the PM-computed eigenvalues*)

GraphicsRow[{ListLogPlot[Abs[eigs/exactish-1 ],PlotRange->All,opts, PlotLabel->"Relative Error vs n"],ListLogPlot[Abs[eigs-exactish],PlotRange->All,opts, PlotLabel->"Absolute Error vs n"]},ImageSize->12*72]

for square well the levels vary roughly quadratically vs n

GraphicsRow[{ListPlot[Abs[eigs-V0],opts],ListLogLogPlot[Abs[eigs-V0],opts]},ImageSize->12*72]