Find eigenvalues and wavefunctions of time independent Schrödinger equation 
for 1D, central, and other separable potentials by The Phase Method

{Now[ ],$Version}

Preamble -- click this opener =>

Grant B. Bunker,  Department of Physics, Illinois Institute of Technology,  bunker@iit.edu    12/25/2024
reference:  “The Phase Method — for numerical solution of the Schrödinger equation”,
 2024, Grant Byrd Bunker, Publisher, Kindle Direct Publishing 
Caveat emptor - this code works well for me, but no warranty is implied. Feel free to write your own using the principles explained in the book. 
Feel free to adapt the code to your own purposes but please cite the book and this code at this URL.

For an initial demonstration using this notebook, in Mathematica (“Wolfram”), just choose Evaluate Notebook from the Evaluation menu. 
Use the Show/Hide Code Cells item from the created palette to switch off/on visual clutter.  Find Evaluating Cell jumps to the cell that is currently executing. 
I recommend that you keep a read-only master template version of this file, and then copy and modify it for different potentials/parameters. 
On Mac OS it’s helpful to use the “stationery” and “locked” attributes for this notebook in the Finder that are accessed from Get Info. 
 You can expand the following cells either by double-clicking their enclosing cell brackets, or by clicking the   “≫”  opener, if present. 
 Email me at bunker@iit.edu with subject line “Phase Method” if you have questions.

Instructions: How to Use This Notebook

Essentially all you just need to enter is your potential function and the upper energy limit to use. This controls the number of eigenvalues to calculate. Keep it to a dozen or so at first to save time.
You can control the accuracy by choosing workingprecision to be twice the desired number of decimal digits of accuracy.  Minimum (fastest) is $MachinePrecision, maximum is >100. 
The ranges are automatically estimated via built-in heuristics, but you can extend or shrink them by changing turning point scale factors (tpsf), or you can override those and use your own values. 
Then just let it rip. 

To evaluate the notebook as a whole, choose “Evaluation/Evaluate Notebook” from Mathematica’s menu. 
The Navigate Cells palette that is created upon running the notebook allows you to hide/show the code or text  to reduce visual clutter. 
Click the buttons while it’s running to jump to the various sections of the calculation to check on progress.
Or simply let it run through and inspect the results at the end.  A bunch of files are automatically created in a new folder named after the potential, unless you turn that off by setting the writefiles flag to False. 

If it emits three quick beeps (or whatever your alert sound is) things went OK.   If it does five slow beeps then something went wrong - inspect the output for explanation.  Usually it’s an input typo or thinko.    Set flag nobeeps=True if you don’t want that. 
To assist long batch jobs it also will attempt to speak the elapsed time at the end to let you know it’s done.
set flag speakelapsedtime=False if you don’t want it to do that.   

At the start of each run,  it’s best to start with a fresh environment.  
I suggest that at the start of each run you first:
• quit the kernel (Evaluation/Quit Kernel/Local)
• delete all output cells (Cell/Delete All Output), and 
• close the Navigate Cells palette(s) that the notebook creates, if it’s leftover from a previous run. A new one will be created.
• you can choose where it appears on your screen by changing the parameters in the WindowMargins function.  It will wrap around if placed off-screen.
• set the flags controlling the overall operation of the notebook (e.g. writefiles, skipwavefunctions, noortho)
• define the potential U and the upper energy limit etest
• set the flags (onesidedpotential, bottomless) describing the potential
• inspect the plots of the potential with etest, showing the turning points
• inspect the Phase Plot to make sure the transitions are sharp and not rounded (i.e. they are “digitized”)
• if they are not, adjust xmin, xmax, or adjust turning point scale factors. This should be handled pretty much automatically
• let it rip (choose the menu item Evaluation/Evaluate Notebook)
• if wavefunctions are calculated, also check the cutoff plots,  orthonormality table, and comparison of eigenvalues with expectation values, for any signs of  irregularities.
• recommend starting with a moderate number of eigenvalues (say 10 or 20) of known potentials at workingprecision=$MachinePrecision and going from there.
• output files are placed into a new directory (folder) that is created in the same location as the notebook itself. 

note about cells:
• to execute a cell, put your cursor in it or select its cell bracket, and hit shift-enter
• you can open an close cells by double-clicking on their cell brackets. You hide a code (input) cell by double-clicking on its output cell. Double-click the enclosing cell to unhide it.  
• the currently evaluating cell is shown as a highlighted cell bracket on the right side of your screen
• to jump to the currently evaluating cell, use the menu item Evaluation/Find Currently Evaluating Cell
• you can hide and show cells using Mathematica’s hierarchical cell brackets shown on the right side of the window
• you can hide/show the code (input) cells using the Navigate Cells palette if you want to de-clutter things.
• you can step through the evaluation one cell at a time by choosing evaluate next cell from the Navigate Cells palette. This can be useful when setting up a job.
• to abort a running evaluation, choose the menu item Evaluation/Abort Evaluation or its keyboard equivalent
• scroll down through the output to examine the output plots, values, and commentary during the calculation. Hiding code and text via the palette may help here. 
• once it seems things are suitably set up, let it rip, and go make a cup of coffee or tea.

Some audio is used for feedback.  Three 100 msec “beeps” indicates that all worked fine. Five 300 msec beeps indicates a fatal error that requires pulling the plug.  Two 200  msec beeps indicates a warning or error with a framed printout that should get attention. 
If speakelapsedtime flag is set, it does that at the end.  Or at least it should.

This notebook is freely shared to demonstrate one practical working implementation of the  Phase Method.
Some of it is pretty kludgey, but feel free to improve on it, or write your own from scratch.
See PM book for theory and rationale and implementation notes.
The code, comments, input, output are interspersed to show how they are related to each other.  
This version has all kinds of bells and whistles, but you want to write your own simplified version, or make a Mathematica Package version, go nuts!   This notebook should easily be adaptable to run as a WolframScript at the command line in batch mode if it were desired. 
See the book for a basic Mathematica version and Python port of the eigenvalue search.

Please cite the PM book if you use this method for producing actual work.  A version of this code is supplementary to the book. 
 
Copyright Grant B. Bunker 2024.   Caveat emptor -- no warranty is implied.

Make navigation palette, show change log

make navigation palette

note: you can place this navigation palette wherever you wish by editing the WindowMargins options below
by default it is placed at the top right of your computer screen but feel free to change that

ClearAll["Global`*"];(*clear all variables in Global context*)
screenwidth=("ScreenArea"/.SystemInformation["Devices","ScreenInformation"])[[1,1,2]];
NotebookClose[navigationpalette];(*cleanup palette(s)*)
nb=SelectedNotebook[];
inputcells=Cells[nb,CellStyle->"Input"];
textcells=Cells[nb,CellStyle->"Text"];mvback:=SelectionMove[nb,Previous,Cell];mvfwd:=SelectionMove[nb,Next,Cell];doit:=SelectionEvaluateCreateCell[nb];
NotebookClose[navigationpalette];(*close palette if is left over from previous invocation*)
navigationpalette=
(CreatePalette[
{Button["Show Code Cells",SetOptions[#,CellOpen->True]&/@inputcells;],
Button["Hide Code Cells",SetOptions[#,CellOpen->False]&/@inputcells;],
Button["Show Text Cells",SetOptions[#,CellOpen->True]&/@textcells;],
Button["Hide Text Cells",SetOptions[#,CellOpen->False]&/@textcells;],Button["Find Evaluating Cell",FrontEndTokenExecute["FindEvaluatingCell"]],
Button["Define Potential etc",NotebookFind[nb,"Define Potential etc",All,CellTags]],
Button["Override energy search range",NotebookFind[nb,"Override energy search range",All,CellTags]],
Button["Override xmin and xmax",NotebookFind[nb,"Override xmin and xmax",All,CellTags]],
Button["homogeneous k>0 est# bound states",NotebookFind[nb,"homogeneous est # bound states",All,CellTags]],
Button["Find Potential Plot",NotebookFind[nb,"Potential Plot",All,CellTags]],
Button["Find Phase Plot",NotebookFind[nb,"Phase Plot",All,CellTags]],
Button["Find Sifting",NotebookFind[nb,"Sifting",All,CellTags]],
Button["Find Refining",NotebookFind[nb,"Refining",All,CellTags]],
Button["Find Eigenvalues",NotebookFind[nb,"Eigenvalues",All,CellTags]],
Button["Find Sanity Checks",NotebookFind[nb,"Sanity Checks",All,CellTags]],
Button["Find Level Plot",NotebookFind[nb,"Level Plot",All,CellTags]],
Button["Find Homogeneous k>0 Check",NotebookFind[nb,"Homogeneous k>0 Check",All,CellTags]],
Button["Find Cutoff Plots",NotebookFind[nb,"Cutoff Plots",All,CellTags]],
Button["Find Orthonormality Table",NotebookFind[nb,"Orthonormality Table",All,CellTags]],
Button["Find Wave Functions",NotebookFind[nb,"Wave Functions",All,CellTags]],
Button["Find Expectation Values",NotebookFind[nb,"Expectation Values",All,CellTags]],
Button["Go To End",NotebookFind[nb,"The End",All,CellTags]],
Button["Move Forward",mvfwd],
Button["Move Back",mvback],
Button["Evaluate Next Cell",{doit,mvfwd}]},
WindowMargins->{{.9*screenwidth,Automatic},{Automatic,0}},WindowTitle->"Navigate Cells"]);mvfwd;

code revision log

==>change Nov 27, 2025
fixed bug that would occur if user mistakenly set the bottomless flag to true for l>0 attractive singular central potentials such as coulomb.  An error like that is caught and the flag reset to false (with message), but other code involved in automatically choosing xmin had to be moved after this reset point to work correctly.  So it was moved.   Also changed override integration method for singular potentials from LSODA to StiffnessSwitching. 
==>changes Oct 25, 2025 
improved wavefunction cutoff code - changed default integrationmethod (ODE solver) f for wave functions from “LSODA”  to “StiffnessSwitching”)
Change default integrationmethod for eigenvalues to StiffnessSwitching which seems to be generally faster and more accurate than LSODA. 
Changed default epsx to 10^(-workingprecision/2) for singular potentials e.g. coulomb, epsx is the minimum x value, at which wavefunction set to zero. 
Change default workingprecision to 60 (to give 30 good digits).
changed ngooddigits to simply Floor[workingprecision/2] (instead of also subtracting 1). 
corrected semiclassical estimates to work correctly when inverse mass scale factor γ  is not 1
=========
changes from Aug 8 2024
	added code revision log 
	added code to verify whether highest wavefunction is actually converged to near zero over the range xmin, xmax, issue warning message if not
	changed default value of tpsfonesided to {10^(-15), 5} 
	adjusted maxbisections (for sifting) and maxiterations (for refining) limits
	added code to prevent double-counting of smallness if nonresidential and bottomless flags are both set
	added workaround for minor graphics error on cutoff plots
	relocated the cell pointed to by “Phase Plot” cell tag so it’s right below actual plot
	added Find Eigenvalues item to Navigate Cells palette with corresponding tag in code
	enhanced warning message for singular potentials U ~ -r^k for k ≤ -2, sea monsters go there
	added warning that epsx affects ultimate accuracy for singular potentials like Coulomb (-1/x) - smaller is better as long as greater than minimum real for specificed working precision
	improved (I think) other text cells
	added test and warning message if abs of differences between eigenvalues is less than 10^(-ngooddigits), i.e. rough accuracy of results
	removed different random fuzzing of intervals - now one fuzz for all intervals
	added 3 beeps at the end to let you know the job is done successfully;  5 slower beeps if run was stopped because of anomalous condition and no point in continuing
	improved  (I think) exception handling and error messages if parallel computation produces duplication of eigenvalues
	for homogeneous potentials, added code to give better estimates of ν for widely varying k - for informational purposes
	added code to not give confusing output for  homogeneous potentials, if the one-sided potential flag is set to True
	changed ‘deleteduplicates’ code to ‘mergeduplicates’ so that it does cluster analysis if needed and take eigs as mean of cluster
	enhanced homogeneous function code to use approximate (OverBar[(Subscript[ν, k])]) expression from book appendix A (especially needed for very large k)
	added code to reset etest if it accidentally set to exceed the continuum limit
	added third-difference printout and more contextual info as to its utility
	added code to disable orthogonality integrals if their number exceeds user-chosen limit
	added internal “sanity check” comparison for potential Abs[x]/2 odd parity eigenvalues
	added  internal “sanity check”  sanity check for finite square well -70, 3 
	reorganized top page layout
	if one-sided potential, modify output file name tag accordingly (nominally same potential can be one-sided or two-sided)
	corrected timing code to use totalstarttime rather than starttime
	added code to speak elapsed time once job complete

Aug 9 fixed typo in text field
	added “homogeneous” flag test for a warning message regarding 1/x^2 potentials
Aug 29 added code to speak elapsed time at end
Sept 1, 2  added better support and documentation for ImplicitRungeKutta 
	    added maxiterationlimit flag allowing (if False, default) refinement to run to completion using FixedPoint
	    changed default workingprecision to $MachinePrecision
	    rounded ngooddigits to integer value to make NumberForm happy
	    added a flag speakelapsedtime  to choose whether to speak elapsed time at end
	    increased precision settings for calculating wave functions 
	    increased default plot range scaling factors
	    moved the location where turning point scale factors tpsf  are set close to where etest is set 
	    changed file name tag to include integration method that is used
	    put in override to enforce LSODA for bottomless potentials
	    added allokbeeps, warningbeeps, and fatalerrorbeeps 
	    added nobeeps flag to shut them off if not desired
	    added Find Sanity Checks to Navigation Palette
	    slightly improved estimated sifting time progress bar code
	    added Find Sifting, Find Refining, Find Sanity Checks to Navigation Palette
	    improved formatting of some sanity checks
Sept 14  corrected preposterous sign error in exponent for default value of short distance turning point scale factor for one-sided potentials
Oct 9  2024  added check so that eminsearch is properly estimated, even if finite epsx gives finite umin and appears to be not bottomless
Oct 20 2024 added a warning about being sure to multiply centrifugal potential by inverse mass scale factor γ, if γ ≠ 1, and l ≠ 0
Nov 29 2024 fixed bug in avoid-redundant-scaling of xmin  when bottomless flag is overridden
Dec 19 2024 changed default MaxRecursion for NIntegrate to 50 for expectation value checks (Automatic not always enough)). 
		    specified a=3 and V0=-9 parameters explicitly for square well “sanity check”. 
		    added more cell tags and navigation items, and various comments
		    changed default fallback for estimated ground state if numeric estimate fails
		    
Dec 24 2024
allow use of “StiffnessSwitching” for eigenvalues, faster and more stable at high workingprecision
default use “LSODA” for wavefunctions
in the Sanity Checks make sure N function is done to at least enough digits
use $MachinePrecision for orthogonality integrals, for speedup
add tags and palette entries to jump to wavefunctions and orthogonality table etc
change wavefunction calculation to use full working precision rather than ngooddigits
where allwfcutplot calculated, change Table to ParallelTable
memoize (cache) values of getwf[ℰ] for fast re-use
change “remove” to “merge” eigenvalues in text description
for cutoffs and uppertps use ParallelTable
calculate getWF at ngooddigits+ 5 
default plotscalefactors->{2,2}
changed warning message for duplicates re suggested precision
changed error message  threshold for orthonormality to 10^(-3)
added print MinMax orthonormality error
Many edits to explanatory text
added code to override selection  of energy search range, and {xmin, xmax} and some checks on values
added code to correctly align eigenvalues for sanity check potentials when energy search range is manually set and eigenvalues don’t start at n=1
relocated setting of epsx (minimum x allowed) to where turning point scale factors are set
changed the spoken alert at the end to use mathematica speechsynthesize instead of system “speak”
reordered items on navigation palette
position palette according to system screen width
changed ParallelTable to Table for orthonormality table (surprisingly much faster!)
Dec 27 tweaked limits and workingprecision of orthonormality integral

=> Define potential and upper energy limit here.
Also set up the calculation parameters and flags if needed, 
then choose Evaluate Notebook from the menu

Set various global flags controlling execution, and a few important parameters

Print["start of execution date/time: ",Now];
totalstarttime=AbsoluteTime[ ];(*keeps track of total execution time for running whole notebook*)
boldface=Style[#,Bold]&;IRK="ImplicitRungeKutta";(*alias - saves typing*)

The following input cell has some flags that set the overall operation of this notebook. 
Change them as you wish.

speakelapsedtime=True;(*whether to speak the total elapsed time when it's complete*)
nobeeps=False; (*disable audio feedback for errors etc*)
maxiterationlimit=True;(*default False: whether to use finite maxiterations in refining process or just let it run to completion*)
writefiles=True; (*flag determines whether or not to write out the plot files*)
nointerrupt=True;(*this flag allows user to choose whether or not to pause for confirmation before proceeding at various places in calculation*)
skipwavefunctions=False; (*this quits the kernel immediately before wavefunctions etc are calculated, to save time in case they are not needed*)
noortho=False;(*if noortho flag is set to true it will disable the time consuming orthogonormality table calculation*)
levellabels=True;(*this is a flag determining whether or not to put numerical eigenvalues on the energy level plot*)
mergeduplicates=False; (*normally False: merges duplicate eigenvalues if they exist - use this only if necessary if you're sure there are no actual near-degeneracies.*)
plotscalefactors={3,3}; (*these are only for plotting, chooses xscale by scaling the turning points - adjust as needed*)

define “beeps” (taps) as audio feedback - these use your default system beep sound
change these as you wish, or use EmitSound if you want to get fancy
three beeps with short intervals between them means all is OK
five beeps with slow intervals between them means something went wrong - check the printed output for more info

If[nobeeps,beep=Null,beep:=Beep[]];
allokbeep:=Do[(beep;Pause[.1]),3];
warningbeep:=Do[(beep;Pause[.2]),2];
fatalerrorbeep:=Do[(beep;Pause[.3]),5];

Define the specific potential function,  some flags describing it,  set the upper energy limit and working precision

Define potential, set flags, and etest, integrationmethod, workingprecision, adjust tpsfs as needed. 
Then let it rip (Evaluate Notebook)

here are a few examples
Morse Potential parameterization as in Dahl and Springborg, JCP 1988:
U0=-70; α=1; x0=3; U=U0(2Exp[-α (x-x0)]-Exp[-2α (x-x0)]); 
harmonic oscillator, unit ℏ ω: U=x^2/2.
Square well test case
U=If[Abs[x]<3,-9,0];*)
Coulomb with 3D centrifugal potential for whatever angular momentum l you want; 
l=2;  U=-1/x + γ l(l+1)/(2 x^2);  γ ≡ (ℏ^2)/m .   Normally set γ=1 in atomic units, but you should set it appropriately for whatever units you use,and the convention used for writing the Schrödinger equation (e.g. whether the 1/2 is included in kinetic energy term).

U=Abs[x]/2;(*V-shape potential, homogeneous, exact solution known*)
Print["U = ",TraditionalForm[U]];
Plot[U,{x,-10,10},PlotRange->All]

γ=1;(*this is an inverse mass parameter in Schrödinger equation*)
workingprecision=60; (* $MachinePrecision is OK for quick results but usually larger values are better*) ;
Print["Working Precision=",workingprecision];
onesidedpotential=False;(*this flag indicates whether or not the domain is restricted to x≥0, e.g. for 3D central potentials where r≥0*)
bottomless=False;(*this flag indicates whether the potential has a minimum or is "bottomless"  e.g. coulomb potential*) 
(*NOTE two-sided bottomless potentials likely would incur divide by zero errors: you can prevent this by clipping the potential to a very large negative value*)
recenter=False;(*determines whether the automatic selection of xmin, xmax using turning point scale factors are to be referenced to the potential minimum rather than x=0*)

User should set this uppermost energy etest here to see what the phase plot (see below) looks like at that energy - this value will be used subsequently, if  {eminsearch, emaxsearch} are not set manually . If a large number of eigenvalues is to be calculated, it may be expedient to set noortho=True  (the number orthogonality integrals goes ~ neigs^2) or even skipwavefunctions=True (if wavefunctions are not needed). Note:  if eminsearch and emaxsearch are manually set by overriding the automatically determined values, etest does not affect the sifting and search processes

etest=4 ; (*this is the user-defined upper energy limit*)

If you want to manually override the automatically estimated energy search range values,  just set the overridesearchrange flag to True,  and enter your own values here, otherwise set the flag to False and let the program calculate it automatically based on the upper energy limit etest and the potential minimum if it exists, and uses heuristics if it doesn’t.

overrideenergysearchrange=False;(*setting this to True disables automatic search range*)
{OReminsearch,ORemaxsearch}={-5,-.001}

Adjust turning point scale factors as necessary. Usually nominal values are OK unless very high accuracy is needed in which case convergence checks vs xmin and xmax should be done. 
They never have to be fine-tuned. TPSFs for potentials like harmonic oscillator or quartic potentials might be {2, 2} or {3, 3}, it doesn’t matter very much as long as they’re several times greater than 1.
The shape of the phase plot will tell you at a glance whether the range is suitable.   You should have sharp transitions, and a region of saturation of the phase plot at the beginning and especially the end.  
If that isn’t the case, extend the x-range. For one-sided potentials you may want to increase the slope boundary condition at xmin by several orders of magnitude. 
Values like {1.2, 1.2} are reasonable for very hard symmetric potentials that keep growing away from the minimum, but don' t have a continuum limit,  e . g . x^12;  
For states near their continuum threshold, such as any barely-bound states of a finite square well, will need larger values of TPSFs, e.g. 10 or 20 or 50 or 100. 
The ultimate accuracy of eigenvalues for singular potentials such as Coulomb is also regulated by the short distance cutoff xmin.  
epsx is set explicitly for bottomless singular potentials - it imposes a lower limit to prevent divide by zero errors. You should change it to check sensitivity of precise eigenvalues.s
Note: for singular potentials, if epsx is too small to be represented by the designated workingprecision, |U[xmin]| effectively may not be numerically representable, and give an error.

epsx=10^(-Round[workingprecision/2]);(*xmin value is explicitly set to epsx for bottomless potentials such as coulomb l=0*)
Print["epsx=",N@epsx];
tpsfonesided={epsx,15};(*nominal turning point scale factors for one-sided potential - adjust as necessary*)
tpsftwosided={5,5};  (*nominal turning point scale factors for two-sided potentials*)

Pick the ODE solver “integrationmethod” here.  LSODA  “old reliable” seems to be the most robust for general potentials, including ones like -1/x. 
But StiffnessSwitching can be considerably faster and seems to work better at high (>40) workingprecision than does LSODA, but is twitchier.
LSODA still seems to be more reliable for Phase Method wave functions once the eigenvalues are accurately known, and that is the default. 
I suggest initially testing with known potentials that are addressed by the internal Sanity Checks in this notebook e.g. Abs[x]/2;  x^2/2; If[Abs[x]<3,-9,0];  -1/x. 
These all have exact (or as-accurate-as-you-want-to-calculate) solutions for verifying results.

integrationmethod="StiffnessSwitching";(*integrationmethod="StiffnessSwitching";*)

do an initial rough plot to check the potential is sensible - in this notebook, xmin and xmax are scaled from the turning points below, although you can override them

{xplotmin,xplotmax}={-20,20};(*temporary values just to take a look - change these as you wish*)
If[onesidedpotential,xplotmin=10^(-3)];(*for plotting only, if restricted to positive x only*)
GraphicsRow[Plot[{U,etest},{x,xplotmin,xplotmax},PlotPoints->200,PlotRange->#,Frame->True]&/@{Automatic,All}]

If[!bottomless &&Minimize[U,x][[1]]==-∞,Print[Framed@Framed@"U apparently has no Minimum - resetting bottomless flag to True"];warningbeep; bottomless=True]//Quiet;

If[!bottomless&&etest<NMinimize[U,x][[1]],(fatalerrorbeep;
ChoiceDialog[Framed@Framed@"Warning: etest/emaxsearch seems to be lower than potential minimum",{"Continue"->"OK","Quit"->"Quit"}])];
If[%=="Quit",Quit[ ]];

(here estimate number of bound states if potential is homogeneous of degree k>0)

Check to see if  potential is homogeneous (this is independent of the Phase Method calculation, but it is useful to know).
If k>0 our semiclassical expression can be used to estimate the number of bound states below a given energy. 
Writing U=(ϵ/2)|(x/a)((|)^k) with ϵ and a energy and length scale parameters, we define η =  (ℏ/( a (Sqrt[m (Abs[ϵ])]))) and μ = ((2 k)/(2+k)).  
Here we are using units in which ℏ->1, m ->1. 
The easily evaluated function  Ω[ℰ_]:= Floor[((2 ℰ/ϵ)^(1/μ) (1/k)!)/((η Sqrt[π] (1/μ)!)) + 1/2]  
i.e.⌊(((((((2 ℰ )/ϵ)))^(1/μ))((1/k))!)/((η (Sqrt[π]))((1/μ))!))+(1/2)⌋
gives the approximate number of bound states below a specified energy ℰ. This can used when planning a calculation for homogeneous potentials with k>0, for example, if you want to set an upper energy cutoff that has a specified number of computed bound states, or to compute a particular range of n, or want to keep the number of states to compute reasonable. 
example computation:  Ωest //. {k -> 1, μ -> 2 k/(2 + k), η→ 1, ϵ→ 1, ℰ→ 10} gives 38 states
The Phase Method calculation of course gives much more accurate results. 

First check to see if potential is homogeneous of degree k.  If it is, writing  U =  (ϵ/2) |(x/a)((|)^k),  
then each eigenvalue then scales exactly as ϵ  ((((ℏ/(a (Sqrt[m|ϵ|])))))^μ)with μ = ((2k)/(2+k)).  
These can be very useful for getting, checking, and estimating eigenvalues; even guesses for shooting method calculations.
set k->2, ν->1/2 for harmonic oscillator,  k-> ∞, ν->0 for infinite square well for which the formula is exact.
ν->1/2 gives same results as JWKB approximation. 
See appendix A of GB book “The Phase Method” for derivation of exact scaling relation and approximate semiclassical formulae.

(*maybek = Simplify[x  D[U, x]/U]*)(*too simple*)

maybek=FullSimplify[Log[(U/.{x->α x})/U]/Log[α],Assumptions-> {α>0}];(*this version handles Abs too*)
homogeneous=NumericQ[ maybek];

the following cell attempts to determine the generic parameter values from the form of the given potential function; 
with Abs etc in the mix,  in some cases it can be trickier than it seems like it should be.

boldface=Style[#,Bold]&;
If[homogeneous && maybek>0,
(Clear[k,μ,η,ℏ,m,a,ϵ,ν,ℰgb];
Print[Framed["Potential "<>ToString[TraditionalForm[U]]<>" appears to be homogeneous of degree k="<>ToString[TraditionalForm[maybek]]]];
ℰgb[n_]:=ϵ/2((((n-ν)η Sqrt[π](1/μ)!/(1/k)!))^μ); (*general form for k>0, ν depends weakly on k; see book*); 
Print[Framed["Approximate expression for energy levels :\n ℰgb[n]="<>ToString[TraditionalForm[ℰgb[n]]]]];
k=maybek; μ=2 k/(2+k);η=ℏ/(a Sqrt[m ϵ]);
savedassumptions=$Assumptions; $Assumptions={a>0,x>0,ϵ>0};
Uform=ϵ/2 Abs[x/a]^k;psoln=Solve[FullSimplify[(0==Log[Uform/U])/.{Abs[s_]->s}],a][[1]];
Print["Potential parameters ","k=",k,"  μ=", μ,"  η=",η];Print["{U, Uform}=",{U,Uform}];(* Print["psoln=",psoln]*);
If[(U==Uform/.psoln//FullSimplify),Print["mapping apparently successful: ",psoln//FullSimplify],Print["apparently failed"]];
     ℰgb[n_]=( ℰgb[n]/.psoln/.{m->1/γ,ℏ->1})//FullSimplify;(*specialized to this case; γ is inverse mass scale factor used in Schrödinger equation*) 
Print[Framed["For this potential, approximate semiclassical energy level estimates are\n(note: JWKB approximation has ν=1/2 (exact for harmonic oscillator);\nsee book for handling much harder and softer potentials):\n
 ℰgb[n]="<>ToString[TraditionalForm[ℰgb[n]//FullSimplify]]]];
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 ℰ/ϵ)^(1/μ)(1/k)!)/((η Sqrt[π](1/μ)!))+νest];
Print[Framed["for this k, but (various ϵ, ℏ, m, a), (Subscript[Ω, est])(ℰ)≈"<>ToString[TraditionalForm[Ωest//.{k->maybek,μ->2 k/(2+k)}]]]];
Ωest=Floor[((2 ℰ/ϵ)^(1/μ)(1/k)!)/((η Sqrt[π](1/μ)!))+νest]/.psoln/.{ℏ->1,m->1/γ}//FullSimplify;
Print[Framed["for this potential with the given parameters (Subscript[Ω, est])(ℰ)≈"<>ToString[TraditionalForm[Ωest//.{k->maybek,μ->2 k/(2+k)}]]]];
Ωest=Ωest//.{ϵ->1,η->1/Sqrt[γ],k->maybek,μ->2 k/(2+k),ℰ->etest};
Print[Framed[boldface["\[Implies] There are an estimated "<>ToString[Ωest]<>" states below etest="<>ToString[etest]]]];
If[Ωest>100,Print[Framed@boldface@"That's a lot of states - you may want to consider a different etest or search range"]];
$Assumptions=savedassumptions;
)]

⇑ Ω estimate ⇑

=> This code  calculates the eigenvalues,  wavefunctions, does various sanity/accuracy checks, and more

If[NumericQ[Ωest],Print[Framed[boldface["\[Implies] Using approximate semiclassical formula (eqn A.3 in PM book): there are an estimated "<>ToString[Ωest]<>" states below etest="<>ToString[etest]]]];]

Define the form of the Schrödinger equation that is used throughout;  If  dimensional units are to be used this is the place to change it. 
This would amount to changing the value -1/2 to  -(ℏ^2)/(2 m) in units natural to the physical system.
Changing the value of the parameter gamma γ above (normally set to 1) will suffice for most purposes.
The resulting coefficient should be on the order of (.1—10) or so if the units are appropriate. 
Don’t use SI units for electron masses etc, the PM can use the dynamic range for its computations.

SetDirectory[NotebookDirectory[ ]];(*this sets the default directory to the location of ths notebook*)
boldface=Style[#,{Bold,FontSize->12}]&;
nf[ndig_Integer]:=NumberForm[#,ndig]&;(*truncate display of digits*)
usewhenevent=True;(*this flag (True normally is the simplest) determines whether ODE solver finds transitions in Phase Plot. Otherwise it's done by other means*)
If[usewhenevent,
Schreq[e_]:={ -γ/2 Φ''[x]+ (U-e)Φ[x]==0,Φ[xmin]==Φ0,Φ'[xmin]==Φp0,WhenEvent[Φ[x]==0,AppendTo[zeros,x]]},Schreq[e_]:={ -γ/2 Φ''[x]+ (U-e)Φ[x]==0,Φ[xmin]==Φ0,Φ'[xmin]==Φp0};](*boundary conditions are numerically set later*)

If[integrationmethod=="ImplicitRungeKutta",ngooddigits=Max[Floor[workingprecision],20],ngooddigits=Round[workingprecision/2]];
(*these estimates are rough and empirical: 
for LSODA and StiffnessSwitching the # of accurate digits reliably approximates (and often slightly exceeds) workingprecision/2 - your mileage may vary for other ODE solvers.  Sufficient data range also is crucial for accuracy*)

Do some checks - analytically (if possible) determine: is it unbounded from below (“is bottomless”)?  
Does it have a continuum limit?  If so, what is its value?

Clear[umin,uminx,uminxsoln];
If[onesidedpotential,xlow=0,xlow=-∞];
tmp=Minimize[{U,x>=xlow},x]//Quiet;(*fallback to see if Minimize fails to find analytical result*)
If[Head[tmp]!=Minimize,{umin,uminxsoln}=tmp;uminx=x/.uminxsoln];
If[umin==-∞,(Print["Potential appears to have no minimum (is bottomless);\nsetting bottomless and onesidedpotential flags to True"];bottomless=True;onesidedpotential=True)];
If[!NumericQ[umin],{umin,uminxsoln}=NMinimize[{U,x>=xlow},x]//Quiet];(*numerical fallback*)
If[onesidedpotential,continuumlimit=Limit[U,x->∞],continuumlimit=Min[Limit[U,x->-∞],Limit[U,x->∞]]];
Print["umin=",umin//N,",  uminx=",uminx];Print["continuumlimit is at ℰ=",continuumlimit];

If[etest>= continuumlimit,(warningbeep;Print[Framed@Framed@"Warning: etest > continuumlimit; resetting etest to (continuumlimit - 0.01)"];etest=continuumlimit-0.01;
Plot[{U,etest},{x,xplotmin,xplotmax},PlotPoints->200,PlotRange->{Max[umin,-10],Min[continuumlimit+(continuumlimit-umin)/3,10]},Frame->True])]

tripwire=(onesidedpotential && overrideenergysearchrange &&((OReminsearch>=continuumlimit) ||(ORemaxsearch >=continuumlimit)));

If[tripwire,(fatalerrorbeep;
ChoiceDialog[Framed@Framed@"Error: energy search range extends into continuum: check eminsearch emaxsearch override",{"Quit"->"Quit"}])];
If[%=="Quit",Quit[ ]];

If[NumberQ[umin] && bottomless,Print[Framed@boldface@"U has a minimum - resetting bottomless flag to False"];bottomless=False;]

If[bottomless &&onesidedpotential,tpsfonesided[[1]]=1];(*prevent redundant rescaling when xmin is explicitly set to epsx*)

If[bottomless &&integrationmethod!="LSODA" && integrationmethod != "StiffnessSwitching",
(warningbeep;
Print[Framed@Framed@"Warning: LSODA or StiffnessSwitching are needed for handling bottomless potentials - resetting integrationmethod to StiffnessSwitching"]; integrationmethod="StiffnessSwitching"; )]

warnlimit=Limit[x^2 (U),x->0,Direction->"FromAbove"];
If[homogeneous&&(warnlimit==∞||(NumberQ[warnlimit]&&warnlimit<0)), (warningbeep;Print[Framed@Framed@boldface@"Warning: attractive potential ∝ -r^k for k≤-2 is full of surprises!\nThere is no lower limit to bound state energies if the potential allowed to extend to x=0!"])]

This does some rough estimates of the length and energy scale of the ground state wavefunction.
It’s most useful for bottomless potentials, but also is used to set the δx grid for catching transitions (if Mathematica WhenEvents are not used).
If there is a potential minimum, find the wavelength near the bottom, then use a grid spacing that is a fraction of that.
If there is no potential minimum, as for coulomb potential, use feynman’s estimate of the scale and energy of the ground state wavefunction.
δx is irrelevant if “usewhenevent” flag is True, because the ODE solver keeps track of zero crossings automatically. This is generally the preferred method. 
the following simple estimate of ground state will fail for x^k  potentials for k≤-2 however

If[bottomless&&NumericQ[umin],eminest=umin,{eminest,x0estsoln}=FindMinimum[1/(2 x^2)+U,{x,Random[ ]}]//Quiet]

If[!NumericQ[eminest],eminest=-10^3];(*fallback if both symbolic and numeric solutions fail; will be improved below*)

If[bottomless,δx=Min[.001,Abs[1/4.*(x/.x0estsoln)]],δx=Min[.001,Abs@N[1/4*2π/Sqrt[2(etest-umin)]]]];

Print["rough ground state energy estimate eminest=",eminest]

If[!usewhenevent,Print["δx=",δx]]

If[bottomless && etest <eminest,Print[Framed@boldface@"Warning: etest may be too low to see any bound states"]]

tpsf is a list of two scaling factors for auto-setting the range {xmin, xmax)}.  These two positive numbers that you specify are used to scale the computed turning points at the energy of interest 
for each state, and thereby to set  {xmin, xmax} for the calculation at that energy.  
For two-sided potentials (e.g. harmonic oscillator) they both should be greater than 1, say 1.1-10 (start with 1.5),  the values depending roughly on the potential.  
If you set these to be too large you may get overflows - if you do, reduce them a bit. 
If the Phase Plot doesn’t saturate (transition become sharp), reduce xmin and increase xmax a bit.
For one-sided potentials (e.g. Coulomb potential with l>0) the lower scale factor should be  1 (if it’s a hard potential) or a little less, say 0.8 (if soft potential), but still greater than zero.
This heuristic is rough but it works pretty well. Alternatively you can instead just enter {xmin, xmax} manually if you wish, but remember  it’s the turning points that run the show.

(*set the turning point scale factors tpsf here:*)

If[onesidedpotential,
(Φ0=0;Φp0=1000;tpsf=tpsfonesided),(*⟸ if one-sided*)
(Φ0=.1;Φp0=0.1;tpsf=tpsftwosided)];(*⟸ if two-sided - adjust as needed*)
SetPrecision[{Φp0,Φ0},workingprecision];   (*this is for the case when extended precision is used*)

if symbolic minimum search etc failed, fall back to numeric search

If[!NumberQ[umin]&& umin != -∞,({umin,uminxsoln}=NMinimize[U,x,Method->"DifferentialEvolution"]//Chop; uminx=x/.uminxsoln;)];  
(*if symbolic gymnastics above fails, try to find global minimum of potential numerically*)
Print["umin=",umin]

Print["etest=",etest,"  umin=",umin//N];
pltminscale=Boole[!onesidedpotential];If[NumericQ[umin],uminplot=umin,uminplot=Automatic];
GraphicsRow[Plot[{U,etest},{x,-20*pltminscale,20},AxesLabel->{"x","U"},PlotRange->#,Frame->True]&/@{{uminplot,Automatic},Automatic}]

Here are two ways to find all of the turning points at a given energy, depending on the nature of the potential. 
There can be many if there are many minima.
If the simple direct NSolve fails, we fall back to a numerical search based on simple mapping and search for sign changes, then linear interpolation over the zero crossing.
 The latter method is a fairly quick and robust way to approximate and map out all the roots of U-ℰ when you can’t solve for them analytically.
 Make sure  extremal range (extmin, extmax) is large enough though.  
 Ultimately we use the extremal values (min, max) of all the turning points. 
Note: the size of the epsx parameter can limit the ultimate accuracy of eigenvalues for bottomless singular potentials - but if it’s too small, computations can get unstable.

If all of this automated rigamarole fails, just choose xmin, xmax manually

timelimit=15;(*try to do accurate numerical search for turning points, but abort if it takes too long*)

Clear[turningpoints0];
turningpoints0[e_]:=turningpoints0[e]=(tplist=x/.TimeConstrained[NSolve[{U==e,x>xlow},x,Reals],timelimit,x->"\[Implies] Abort: accurate turning point search exceeding time limit "<>ToString[timelimit]<>"sec"]//Quiet;If[onesidedpotential,tplist=Select[tplist,#>0&]];(*take only positive tps*);
If [onesidedpotential && (Length[tplist]<2),tplist=Join[{epsx},tplist]];
If [!onesidedpotential && (Length[tplist]<2),tplist=Flatten@Join[-tplist,tplist]];
Return[tplist]
);(*this is the simplest method if it works*)

Off[General::munfl];(*shutoff underflow warning message*)

(* otherwise use this is a brute force method (this seems good for potentials with many minima, but have to ensure extremal range (-1000,1000} is large enough)*)
{extrxmin,extrxmax,dx}={Max[xlow+epsx,-1000],1000,.1};(*limits for x-ranges over which to search for turning points if using brute force method*)
Clear[Ufun]; Ufun[x_]=U; 
fuzz=10^(-ngooddigits)*RandomReal[{.1,1}];(*to try to avoid cases where an exact eigenvalue happens to be at a energy interval endpoint*)
{xray,uray}=Transpose@Table[{x,U},{x,extrxmin,extrxmax,dx+fuzz}];xray=Drop[xray,1];
rootintervals[e_]:={xray[[#-1]],xray[[#+1]]}&/@Flatten@Position[Abs@Differences[Sign[uray-e]],__?Positive];
root[{{x1_,y1_},{x2_,y2_}}]:=x1-y1((x2-x1)/(y2-y1));(*simple slinear interpolation for approximate root*)
turningpoints1[e_]:=(rooti=rootintervals[e]; intrvl=(Transpose[{#,e-Ufun/@#}])&/@rooti; tplist=root/@intrvl; 
If [bottomless&&Length[tplist]<2,Return@Join[{epsx},tplist]];
Return[tplist])

(*use method 0 if it works, otherwise use method 1*)
turningpoints=.;tps=turningpoints0[etest];
If[Length[tps]>=2&&(And@@NumberQ/@MinMax[tps]),turningpoints=turningpoints0,turningpoints=turningpoints1];

Clear[extremalturningpoints];
extremalturningpoints[e_]:=turningpoints[e]

we want to make sure xmin is at a low enough value that the wave function has “time” to saturate for small x
(jacking up the initial slope and value can get it to saturation faster)
also make sure xmax is large enough that the wave function has “time” to “collect” all the eigenvalues with newly appearing transitions that first appear at large x
we will make an automated estimate of a good xmin, xmax by finding turning points at highest energy to be calculated, and going a bit more than that away from the minimum
the Phase Plot will tell you when you’ve gone far enough; 
overflows and underflows will tell you if you have gone too far.
So here we will specify the integration range by extending the classical turning points by a scale factor

Print["{onesidedpotential, bottomless}=",{onesidedpotential,bottomless}]

Print["{etest, umin}=",{etest,umin}]

tplist=Sort[MinMax[turningpoints[etest]]];

Print["turning point list: tplist=",tplist]

(* If all else fails*)

tpfallback={-30,30};
If[!NumericQ[tplist[[1]]],tplist[[1]]=tpfallback[[1]];Print["Warning: Resetting non-numeric turning point 1 to nominal value "<> ToString[tpfallback[[1]]]]];
If[!NumericQ[tplist[[2]]],tplist[[2]]=tpfallback[[2]];Print["Warning: Resetting non-numeric turning point 2 to nominal value "<> ToString[tpfallback[[2]]]]];
If[onesidedpotential &&bottomless,tplist[[1]]=epsx];(*override for one-sided potential*);

Print["{U, umin, etest, tplist}=",{U,umin,etest,tplist//N}];

Print["turning point scale factors=",tpsf];

If[onesidedpotential ||(!recenter),{xmin,xmax}=tpsf*tplist,
{xmin,xmax}=tplist+tpsf*(tplist-Mean[tplist])];(*handle offset minimum in potential if present*)

If[onesidedpotential, xmin=Max[xmin,epsx]];

if you want to manually override xmin and xmax,  set those values here
 probably better to use scaled turning points though

(*{xmin,xmax}={10^(-15),10^3}*)

set precision of input variables so they don't spoil the precision of the result

{xmin,xmax}=SetPrecision[{xmin,xmax},workingprecision];
Print["{xmin,xmax}=",N@{xmin,xmax}]

the dots on the plot below show the turning points for energy etest

Print["umin, etest, xmin, xmax=",{umin,etest,xmin,xmax}]

tpxbar=Mean[tplist];(*used for offset minimum in potential*)
tpdeltax=tplist[[-1]]-tplist[[1]];
If[onesidedpotential,tplist[[1]]=epsx];

PotentialPlot

If[onesidedpotential,{xplotmin,xplotmax}=plotscalefactors*tplist,{xplotmin,xplotmax}=tpxbar+plotscalefactors[[2]]*(tplist-tpxbar)];
Print["{xplotmin,xplotmax}=",{xplotmin,xplotmax}]
potentialplot=Plot[{U,etest},{x,xplotmin,xplotmax},Frame->True,PlotRange->All,Epilog->{PointSize[.02],Point[Transpose[{tplist,{etest,etest}}]]}]

Print["onesidedpotential=",onesidedpotential]

Phase Plot

Check the Phase Plot below, evaluated near the upper energy where the eigenvalues are sought 
if this looks good (“digitized”, you’re good to go - the rest can probably run on its own)

Here is a Phase Plot to check on the range for the given energy. 
It’s a good idea to do a check at the upper end of the energy range you want to compute.
we want the Phase Plot to saturate promptly and to have sharp steps in between, and an extended saturation region at the end.
If features are rounded or simply oscillatory, lower xmin and raise xmax so as to get saturation and step,  or slightly reduce the upper energy limit to be calculated.
If overflow or underflows bother you,  increase xmin and/or lower xmax, but they are usually not a problem, unless it is excessive so the solution constantly is hitting the ceiling. This can occur for potentials that grow rapidly vs x to very large values when setting turning point scale factors to nominal values e.g. {2,2}.  In such cases the cutoffs and calculated wave functions likely be visibly incorrect (eigenvalues, interestingly, not so much).   Orthonormality tables should tell you of the problem if it occurs. In such cases the turning point scale factors simply will need to be reduced to values closer to 1, e.g. {1.1, 1.1}, and then it should run fine.

If[!usewhenevent,Print["δx=",δx]](*est minimum sampling needed to resolve transitions*)

Clear[Φ];zeros={};
Print["umin=",umin,"  energy=",etest];
Print[AbsoluteTiming[(Off[General::munfl,General::unfl,General::ovfl,InterpolatingFunction::dmval];soln1=(NDSolve[Schreq[etest],Φ[x],{x,xmin,xmax},Method->integrationmethod,WorkingPrecision->workingprecision,MaxSteps->10^7]//Quiet)[[1]];)][[1]]," sec"];
phaseplot[x_]=Tanh[Φ[x]/.soln1]//Quiet;
If[usewhenevent,nboundstates=Length[zeros],
(transitions=DeleteCases[0]@Differences@Table[Sign[phaseplot[x]],{x,xmin,xmax,δx}]/2;nboundstates=Length[transitions];)];
phaseplt=Plot[phaseplot[x],{x,xmin,xmax},PlotRange->All,PlotPoints->2000,Frame->True,FrameLabel->{{"tanh(Φ(x))",None},{"x",boldface@("Phase Plot: energy="<>ToString[TraditionalForm@etest]<>";  Estimated # bound states= "<>ToString[nboundstates])}}]

If[!(NumericQ[nboundstates]&& nboundstates>0),(fatalerrorbeep;
ChoiceDialog[Framed@Framed@"Error: for some reason, no bound states were found below etest - you need to quit now and figure out the problem",{"Quit"->"Quit"}])];
If[%=="Quit",Quit[ ]];

⇑  Phase Plot  ⇑

the following interrupt is to allow you to inspect the Phase Plot above for reasonableness
it will execute if the flag nointerrupt is set to False
if transitions aren’t sharp, click Abort, and extend the range, and rerun it
if it’s OK click Continue

If[!nointerrupt,(NotebookFind[nb,"Phase Plot",All,CellTags];Interrupt[ ])]

If you want to see it, the following line of code generates a list of transitions (sign changes) that occur in the Phase Plot  above
(the length of which is the number of zero crossings which estimates the number of bound states).
For this purpose, we don’t care about the locations in x of the transitions, only their existence and number.

(*transitions=DeleteCases[0]@Differences@Table[Sign[phaseplot[x]],{x,xmin,xmax,.01}]/2;*)

Print["{U, umin, etest, tplist}=",{U,umin,etest,tplist//N}]

Define functions to find eigenvalues

Look for eigenvalues first by sifting, then refining

the nbs function (“number of bound states”; referred to as Ω(ℰ) in book) returns two values: the energy and the # of bound states below that energy

following is a binary search  that bisects the given interval,  and generates new subintervals according to whether or not it finds a change in the number of bound states between the endpoints.  
Each bisection takes an input interval, and produces either the old interval or 1 or 2 new intervals, depending on whether it found zero or  nonzero number of eigenvalues. 
Any interval with no eigenvalues in it is discarded.  Each iteration of bsearch should only increase the number of intervals by one or two at most, until all the eigenvalues within the original interval are acquired.    
After all intervals containing a single eigenvalue are determined (“sifting”), their further “refinement” is done in parallel. 
Refinement is repeated until a fixed point is found, or (as a hedge) a maximum number of iterations is exceeded. That  maximum is chosen to give the desired number of good digits (typically half of the working precision).

Print["{umin, etest, xmin, xmax, workingprecision}=",{umin,etest,xmin,xmax,workingprecision}]

nbs determines the number of bound states below the specified energy by counting the transitions in the Phase Plot.
Be sure to execute the following ClearAll[nbs] if you want to clear out the cached (memo-ized) previously computed values

If[NumericQ[umin],Print["(Subscript[λ, min])=",2π/Sqrt[2(etest-umin)]//N]](*this is a check to make sure sampling in x is fine enough to catch transitions *)

nboundstates=.;
ClearAll[nbs];(*clears out previously memo-ized values*)
nbs[e_]:=nbs[e]=(zeros={};soln=Quiet@NDSolve[Schreq[e],Φ[x],{x,xmin,xmax},Method->integrationmethod,MaxSteps->10^7,WorkingPrecision->workingprecision][[1]];
phaseplot[x_]=Tanh[Φ[x]/.soln];
If[usewhenevent,nboundstates=Length[zeros],
(transitions=DeleteCases[0]@Differences@Table[Sign[phaseplot[x]],{x,xmin,xmax,δx}]/2;nboundstates=Length[transitions];)];
Return[{e,nboundstates}])

this binary search uses the change in the number of bound states between two energies as criterion for eigenvalue in that interval
it bisects the given interval and then returns zero, one or two new intervals that contain one or more eigenvalues. Any interval with no eigenvalues in it is discarded.

Clear[nstates];
nstates[e_]:=nbs[e][[2]];
ClearAll[bsearch];
intervals={};
bsearch[intrvl_]:=
Module[{slow,shigh,smid,mid,newintervals},
({low,high}=intrvl;
mid=Mean[{low,high}];
slow=nstates[low];
shigh=nstates[high];
smid=nstates[mid];
newintervals={ };
If[slow<smid,AppendTo[newintervals,{low,mid}]];
If[smid<shigh,AppendTo[newintervals,{mid,high}]];
Return[newintervals]
);]

Choose the range {eminsearch, emaxsearch} within which to search for all the eigenvalues, and then actually do it

Repeated application of bsearch function eventually generates a list of intervals, each containing exactly one eigenvalue
For eminsearch, here we use umin (the minimum of U, if there is one).
If there is no minimum to U (it’s “bottomless”), compute a rough estimate a la Feynman as described in lectures.
If that fails utterly (e.g. 1/x^2 potential) just put in some a number.  
Eminsearch and emaxsearch can be overridden manually as needed, just change the numbers.

If[!NumericQ[eminest],eminest=-10^3];(*if numeric estimate failed*)

if eminsearch still seems excessively low, do a progressive scan downward in energy, using the bound state counter nbs,
 to a find an energy at which there are no lower energy states.  
 If this quasi-automated rigamarole isn’t satisfactory despite our best intentions
 (anticipating all possible potential forms and variants is a bit complicated),  
 just manually put your own numbers  in {eminsearch, emaxsearch)  where they are first used for to set startinginterval below

escanfactor=10;(*increases by an order of magnitude each time*)

TimeConstrained[
If[eminest<=-10^6,(i=1;et=-.01;(* this is the(negative) energy from which to start scanning downward for determining eminsearch*)
If[et<0,Print["searching for lower limit of states"];While[(nbs[et][[2]]>0 && i<40),Print@nbs[et];
et=escanfactor*et//N; i++]];
eminest=et;Print["# test energies=",i];
Print["eminest=",eminest])],20,(Print["Taking too Long: Abort"];eminest=et;)]

Print["{umin,etest,xmin,xmax,eminest}=",{umin,etest,xmin,xmax,eminest}];

If[maybek<=-2,Print[Framed@"Warning: U ∝ (x^(-k)) for k≤-2 (and l=0)and xmin=0 there is no lower limit to bound state energy; you may need to increase xmin to limit # of bound states"]];
If[bottomless,Print["eminest: rough estimate for minimum energy eigenvalue = ",eminest]];
eminsearch=If[eminest<0&&(umin<-10^3 ),2*eminest,umin];(*factor of two safety factor; this value may be overridden later*)
eminsearch=Max[eminsearch,-10^6];(*clip it if necessary*)
emaxsearch=etest;(* set this upper limit to whatever you want*)
Print["{umin, etest, eminsearch,emaxsearch}=",{umin,etest,eminsearch,emaxsearch}//N]

If overrideenergysearchrange flag is set,  the values set above are used

If[overrideenergysearchrange,{eminsearch,emaxsearch}=Sort@{OReminsearch,ORemaxsearch}];

{eminsearch,emaxsearch}

startinginterval=SetPrecision[{eminsearch,emaxsearch}//Chop,workingprecision];
Print["eminsearch=",eminsearch,",  emaxsearch=",emaxsearch]

timelimit=60;(*secs*)
TimeConstrained[(statesbelow=nbs[eminsearch][[2]];(*this might take a little time depending on potential and etest*)
If[statesbelow>0,Print[Framed@boldface@("There appear to be "<>ToString[statesbelow]<>" states below eminsearch="<>ToString[eminsearch])]];),timelimit,Print["Aborted: evaluating nbs at eminsearch is taking too long"]];

“sift” uses binary search to identify intervals in which there are single eigenvalues 
this generates a bunch of intervals that can be separately “refined” in parallel

define sift and refine functions; using maxiterations is a failsafe against it taking forever to get to fixed point

Off[General::ovfl,General::unfl];

maxiterations=Ceiling[(ngooddigits*Log[2,10])]+10;(*for refining*)
If[!maxiterationlimit,maxiterations=∞];
sift=bsearch/@(Partition[Flatten[#],2])&;
refine=FixedPoint[sift,#,maxiterations]&;

{maxiterations,2.^(-maxiterations),10^(-ngooddigits)//N,10^(-workingprecision)//N}

this does the sifting (recursively divides maxbisections times to separate out the eigenvalues, mostly for visualization, 
but the effort is not wasted for later refinement because previously computed values are cached (memo-ized).
Here this sifting process is done serially, but the later refinement to higher precision is done in parallel on whatever cores are available.
Note that, even if levels are not fully resolved in the sifting process, they should be resolved during later refinement.
If you want to limit the time spent sifting just reduce maxbisections a smaller number say 7.

{thetimings,nbsvals}=Transpose[AbsoluteTiming[nbs[#]]&/@{eminsearch,Mean[{eminsearch,emaxsearch}],emaxsearch}];
Print["eminsearch, emaxsearch=",{eminsearch,emaxsearch}];
Print["Ω[eminsearch], Ω[emid], Ω[emaxsearch]",nbsvals[[All,2]]];

this gets the number of eigenvalues in the search interval and roughly estimates the number of bisections needed to get them all (then pads it a bit)
the following code simply tries to estimate how long it should take

timepereval=Max[thetimings];
maxbisections=Ceiling@Log[2,(nbs[emaxsearch][[2]]-nbs[eminsearch][[2]])]+7;(*for sifting - so as to roughly separate eigenvalues*)
sparsenessfactor=30;(*your mileage may vary*)
(*If[NumericQ[maybek]&& maybek<0,sparsenessfactor=450];(*levels for 1/(r^k) potentials are sparser - affects time estimate*)*)
est=siftingtimeestimate=2^maxbisections/(sparsenessfactor)*timepereval;
endest=est;
timepereval sec

Sifting this does a rough separation of the eigenvalues (“sifting”) which typically takes a few minutes to do as a serial process
(yes, this could be parallelized also)

Print[Framed@boldface@("A very rough predicted (serial) sifting time under these conditions is ~ "<>ToString[NumberForm[siftingtimeestimate/60.,{5,2}]]<>" minutes")];

ProgressIndicator[Dynamic@Clock[est,est,1],{0,Dynamic[endest]}]

starttime=AbsoluteTime[ ];
{sifttime,intervals}=AbsoluteTiming[Flatten[Nest[sift,startinginterval,maxbisections],1]];
endest=0;Print["sifting complete: sifting time=",sifttime," sec,  # of sifted intervals=",Length[intervals]];
ListPlot[Tooltip[Interval/@intervals],PlotRange->All,Frame->True,PlotMarkers->"OpenMarkers",PlotLabel->"`Sifted' rough eigenvalue estimates",FrameLabel->{"n","ℰ"}]

the following figure shows the evolution of the previously completed “sifting” process
it labels the intervals with the number of eigenvalues they contain (yes, it’s kind of a kludge; adjust code as needed)

Clear[spi,thelists,rougheigs,dens,demeans,deneigs]
rougheigs=Mean/@intervals;
nits=Min[maxbisections,7];(*levels of iteration to show in the sifting plot*)

spi=siftplotintervals=Interval/@Partition[Flatten[NestList[sift,startinginterval,nits]],2];
xlabelshift=.015*Length[spi]; ylabelshift=0;(*just a guess: use whatever works well for sifting plot*)
thelists=Flatten[List@@@spi,1];
neigsin[intrvl_]:=Length@Select[rougheigs,IntervalMemberQ[intrvl,#]&];
dens=xlabelshift+Range[1,Length[spi]];(*the n values that set the horizontal position of labels for each interval*)
demeans=ylabelshift+Mean/@thelists;(*the means that set the vertical position of labels for each interval*)
deneigs=ToString/@(neigsin/@spi);(*the number of eigenvalues contained in each interval, converted to strings*)
bf=Style[#,{Bold,Larger}]&;
deepilog=Text@@@(Transpose[{ToString/@deneigs,Transpose[{dens,demeans}]}]);
siftplot=ListPlot[siftplotintervals,PlotRange->All,PlotStyle->{PointSize[.001]},IntervalMarkers->"Fences",IntervalMarkersStyle->{Thickness[.0025]},ImageSize->Large,Frame->True,Ticks->Full,GridLines->All,Epilog->deepilog,FrameLabel->{{bf[(Subscript[ℰ, n])],None},{bf@"sifting sequence",bf["Convergence of search intervals during early sifting process;  U="<>ToString@TraditionalForm[U]]}}]

the number of eigenvalues in the search interval {eminsearch, emaxsearch} that ultimately should be found (after full refinement, at least) is:

eps=10^(-ngooddigits);(*we have to evaluate the number of bound states nbs (AKA Omega) just slightly above an eigenvalue for it to show up*)
AbsoluteTiming[Print["number of states found=",nbs[emaxsearch+eps][[2]]-nbs[eminsearch+eps][[2]]]]

Refining Now we refine the eigenvalues by zooming in on the intervals that were obtained by sifting
Mathematica will present a progress bar and time estimate below for the parallel computation

nomas=Off[General::unfl,General::ovfl,InterpolatingFunction::dmval,LinkObject::linkd,KernelObject::rdead];(*turn off overflow underflow etc  messages - those conditions don't much matter here; and similarly for parallel kernel communication messages*)
ParallelEvaluate[nomas];(*and also do it for parallel kernels*)
separator=Style[StringJoin@ConstantArray["*-",50],Bold];
{refinetime,result}=AbsoluteTiming[ParallelMap[Quiet@refine,intervals]];
Print[separator];
Print["refining complete: refining time=",refinetime,"sec;   # of intervals=",Length[intervals]];
intervals=Flatten[result,2];
elapsedtime=AbsoluteTime[ ]-starttime;
Print["elapsed time since start=",elapsedtime," secs"];
Print["(sifting+refining time)",(sifttime+refinetime) secs];
Print["(sifting+refining time)/eig=",(sifttime+refinetime)/Length[intervals] (secs/eig)];
Print[separator];
ListPlot[Interval/@intervals,PlotRange->All,Frame->True,PlotLabel->"Eigenvalue Intervals from Parallel Refinement",FrameLabel->{"n","ℰ"},IntervalMarkers->Automatic]

Print["mergeduplicates=",mergeduplicates];

nboundstates=nstates[emaxsearch]-nstates[eminsearch]

If[Length[intervals]>nboundstates,Print[Framed@"Warning: refining process may have picked up duplicate eigenvalues;\n consider changing workingprecision"]]

check to verify that we got all the eigenvalues in the range note that the range does not have to start with the lowest eigenvalues
the following code evaluates nbs  (number of bound states aka integrated density of states Ω) at an energy slightly higher than the eigenvalue so as to catch the bound states,  then packages up that state number with the eigenvalue itself

eps=10^(-ngooddigits);

if spurious duplicate eigenvalues are found in the (parallel) refine process, and if the mergeduplicates flag is set to True, the code below will merge the ones that cluster together so as to give the known number of states.  It does this by doing a cluster analysis dividing into exactly nboundstates clusters, where nboundstates is their number the phase plot calculates are between eminsearch and emaxsearch.  This duplication is unusual, but can occur as a consequence of the parallel execution and numeric instabilities. In such cases changing the number of digits of extended precision may suffice to eliminate the problem.
Normally it is best  to leave the mergeduplicates flag set to False in case there are some unexpected degeneracies, and only invoke it if absolutely necessary.

clusters=SortBy[FindClusters[Mean/@intervals,nboundstates],Mean];
Print["estimated raw eigenvalues (or clusters)=\n"];
Print[TableForm@clusters];

if there are any duplicates, the standard deviations of the clusters will be found here:

If[Length[#]>1,DeleteCases[({Mean[#],StandardDeviation[#]}&/@clusters),Null]]

if mergeduplicates flag is set, any such variations will be averaged (not standard procedure).

If[mergeduplicates,eigs=Flatten[Mean/@clusters],eigs=Mean/@intervals];
neigs=Length[eigs];
Print["number of good digits=",ngooddigits];

if user doesn’t want to merge “duplicates”, and there are still discrepancies, shut it down with an error message

If[nboundstates <Length[eigs],(Print[Framed@Framed@boldface@"Error: # of found eigenvalues is greater than the number of bound states from Phase Plot - there may be duplicates: consider changing workingprecision (suggest 25...31) or setting mergeduplicates flag to True (if standarddeviations small enough) ; Quitting Now"];fatalerrorbeep; Quit[ ])];(*5 slow beeps means it didn't go so swell*)

neigs=Length[eigs]

Print["sorted=",eigs==Sort[eigs]](*do the eigenvalues increase uniformly with n as they ought to?*);

eigenvaluelist=Transpose[{eigs,(nbs/@(eigs+eps))[[All,2]]}];
bslist=eigenvaluelist[[All,2]];
nofeigs=Range[bslist[[1]],bslist[[1]]+Length@bslist-1];
(*Print["found eigenvalues"<>ToString[bslist]];*)
If[bslist==nofeigs,
Print[Framed@boldface@"OK - all eigenvalues apparently found"],(*Print[Framed@boldface["found eigenvalues #"<>ToString[bslist]<>";\n apparently missing eigs# "<>ToString@Complement[nofeigs,bslist]<>"\n; (misnumbering is not necessarily a problem)"]]*)];

integrationmethod

(*Print[boldface@("eigenvalues rounded to ngooddigits = " <> ToString[ngooddigits] <> "\n"), Framed@Framed@TableForm[Transpose[{nofeigs, nf[ngooddigits] /@ eigs}]]]*)

Print[boldface@("eigenvalues rounded to ngooddigits = "<>ToString[ngooddigits]<>"\n"),Framed@Framed@TableForm[Transpose[{nofeigs,NumberForm[#,{ngooddigits,ngooddigits},DigitBlock->3]&/@eigs}]]]

eigs==Sort[eigs]

Iconize[bslist](*these are the corresponding state numbers*)

Iconize[eigs](*copy this object into another notebook to use the eigenvalues there, or use the output files*)

Sanity Checks for some specific potentials   =====================================================================================================================

If[integrationmethod!="LSODA"&&integrationmethod!="StiffnessSwitching",Print[Framed@("note: Caveat Emptor regarding ngooddigits for integration method "<>integrationmethod)]];

Abs[x]/2 linear V-shape potential:  calculate relative error of computed eigenvalues compared to exact eigenvalues for odd parity levels

(*bslist is list of bound states numbered from ground state*)
oddparitystatenumbers=Select[(bslist)/2,IntegerQ];
nstart=If[EvenQ[bslist[[1]]],1,2];

If[ U===Abs[x]/2 && !onesidedpotential && γ==1,
(Print["workingprecision=",workingprecision," ngooddigits=",ngooddigits];oddparityeigs=eigs[[nstart;;;;2]]; exactoddparityeigs=-1/2 AiryAiZero[Select[bslist,EvenQ]/2];  relerr=oddparityeigs/exactoddparityeigs-1;  Print[bf@"\nRelative error with respect to exact values for U=Abs[x]/2\n(odd parity levels only, truncated to estimated ngooddigits):\n",nf[ngooddigits]@TableForm[Transpose[{oddparityeigs,N[exactoddparityeigs,ngooddigits],relerr}],TableHeadings->{None,{"Phase Method","  Exact","relative error"}}]])];

x^2/2 harmonic oscillator potential only:  calculate relative error of computed eigenvalues compared to exact eigenvalues 
use of inverse mass scale factor complicates slightly

n0=Round[(eigs[[1]]/Sqrt[γ]+1/2)] (*find which n the lowest eigenvalue corresponds to if doesn't start at 1; gamma is inverse mass scale factor in Schrödinger equation*);

If[ U===x^2/2 && !onesidedpotential,Print[{"WP="<>ToString[workingprecision],"ngooddigits="<>ToString[ngooddigits]}];Print["\nRelative error compared to exact values for U=x^2/2 and γ="<>ToString[γ]<>":"];Print[Framed@TableForm[eigs/Table[(n-1/2)*Sqrt[γ],{n,n0,n0+Length[eigs]-1}]-1]]];

-1/x coulomb potential only:  calculate relative error of computed eigenvalues compared to exact eigenvalues

{nbeg,nend}=Round[Sqrt[-1/(2#)]]&/@{eigs[[1]],eigs[[-1]]}  (*find which n the first and last computed eigenvalue correspond to, in case they don't start at n=1*);

{workingprecision,ngooddigits}

If[U===-1/x && onesidedpotential &&γ==1,(relerrs=eigs/Table[-1/(2 n^2),{n,nbeg,nend}]-1;
Framed@TableForm[Transpose[{Range[nbeg,nend],relerrs}],TableHeadings->{None,boldface/@{"n","Error Relative to Exact"}}])]

this is another sanity check,  for U===If[Abs[x]<3,-9,0]   (square well, depth -9, half-width 3)
compared to iconized rootfinding eigenvalues that were calculated with working precision 100 digits RFeigs=(*100 digit WP*);
For the highest states near the continuum be sure the turning point scale factors are large enough (e.g. >10 or more) for your desired accuracy

If[U===If[Abs[x]<3,-9,0]&&γ==1,({nbeg,nend}=(Position[N[,6],N[#,6]][[1,1]]&)/@{eigs[[1]],eigs[[-1]]};Print[{workingprecision,ngooddigits}];
selectedeigs=Take[,nbeg;;nend];Print[bf@("\nRelative error with respect to precise rootfinding eigenvalues for square well: U(x)="<>ToString[U]<>"\n"),Framed@NumberForm[TableForm@Transpose[{eigs,selectedeigs,eigs/selectedeigs-1}],ngooddigits]])]

End of sanity checks for some specific potentials   =================================================================================================================

Here are first,  second, and third differences.  These can give useful information on the n-dependence of the eigenvalues. 
The first differences are just the level spacings. If the second differences are zero, the levels vary linearly with n (as in harmonic oscillators). 
If third differences are zero, the levels vary  quadratically with n (as for Morse potentials).

first differences  (if these are nearly constant the eigenvalues are approximately (or exactly) a linear function of n)

Print[boldface@"first differences between eigenvalues=\n",nf[ngooddigits]@TableForm[Differences[eigs]]];

second differences (if these are nearly constant the eigenvalues are approximately (or exactly) a quadratic function of n)

Print[boldface@"second differences between eigenvalues=\n",nf[ngooddigits]@TableForm[Differences[eigs,2]]];

third differences  (if these are nearly constant the eigenvalues are approximately (or exactly) a cubic function of n)

Print[boldface@"third differences between eigenvalues=\n",nf[ngooddigits]@TableForm[Differences[eigs,3]]];

here the raw levels are graphed

nofeigs=Range[1,Length[eigs]];(*renumber from the lowest state as needed; sometimes they may be offset by 1*)

Print["state numbering:=",nofeigs]

frmlbl={{"(Subscript[ℰ, n])",None},{"n","Eigenvalue spectrum:  U = "<>ToString[TraditionalForm[U]]}};
loglogfrmlbl={{"|(Subscript[ℰ, n])|",None},{"n","Log-Log |(Subscript[ℰ, n])|:  U = "<>ToString[TraditionalForm[U]]}};
eigplot=Framed@GraphicsRow[{ListPlot[Transpose[{nofeigs,eigs}],PlotRange->All,Joined->True,PlotMarkers->"OpenMarkers",FrameTicks->All,GridLines->Full,Frame->True,FrameLabel->frmlbl],ListLogLogPlot[Transpose[{nofeigs,Abs[eigs]}],Joined->True,PlotMarkers->"OpenMarkers",FrameTicks->All,GridLines->Full,Frame->True,FrameLabel->loglogfrmlbl]}]

this gives a summary plot of the energy levels;  basically it displaces the level eigenvalue labels horizontally if they get too close
(yes, the label position shift stuff is a hack);    set levellabels=False to turn off the labels

levellabels

tps=tplist;

plotscalefactors={1,1.2};(*change these to whatever you like*)
{xplotmin,xplotmax}=plotscalefactors*tps

If[NumericQ[umin],emin=Max[eminsearch,umin],emin=eigs[[1]]];
ediffs=Differences[eigs];ediffs=Append[ediffs,ediffs[[-1]]];(*find the energy differences, pad by last element so same length*)
xs=If[#<(emaxsearch-eminsearch)/70,.15*(xplotmax-xplotmin),0]&/@ediffs; (*if spacing too small between levels, displace label in x *)
Table[lblxshift[eigs[[i]]]=xs[[i]],{i,1,Length[eigs]}];(*make a lookup table of shifts indexed by eigenvalue *);
bars={{xplotmin,#},{xplotmax,#}}&;
lblyshift=(emaxsearch-eminsearch)/100.;(*vertical displacement of level label in plot, choose as you will*)
lbls=Text[ToString[nf[ngooddigits][#]],{xplotmax*.6+lblxshift[#],#+lblyshift}]&;
dele=(eigs[[-1]]-eigs[[1]]);(*for scaling plot*); eminplot=emin -0.05dele ; emaxplot=eigs[[-1]]+.05dele; plotrange={eminplot,emaxplot};
levelplot=Plot[U,{x,xplotmin,xplotmax},PlotRange->plotrange,FrameLabel->{{ℰ,None},{"x",boldface["U="<>ToString[TraditionalForm@U]]}},Frame->True,ImageSize->10*72,Epilog->{Line[bars/@eigs],If[levellabels,Opacity[1],Opacity[0]],lbls/@eigs}]

here are the energy eigenvalues:  the digits should be correct except maybe the last one - increase working precision if you need a few more. 
If you need ultimate accuracy and the potential is smooth you can use Lagrange Mesh method.

neigs=Length[eigs];Print["number of eigenvalues in range = ",neigs]

Print["estimated number of good digits (caveat emptor)=",ngooddigits]

Print["eigs=",nf[ngooddigits]/@eigs]

Iconize[eigs](*here we save eigenvalues as an Iconized object - copy it into another notebook to use the eigenvalues there, otherwise use output files*)

Homogeneous k>0 Check results vs approximate homogeneous potential formula
This compares the calculated PM eigenvalues with our approximate semiclassical formula given A.2 in the Phase Method book appendix A.
Exact values are ν= 0 for infinite square well (which corresponds to k-> ∞);  ν=(1/2) for harmonic oscillators (k= 2). WKB value is ν->1/2. 
Empirical approximate formula for ν from book is used here. It deviates from conventional semiclassical values.  For more accurate estimates at low quantum number n,  and deviations from k=2,  ⋁> (1/2) for 0<k<2, choose ⋁< 1/2 for k>2. 
Note: for these potentials of form U=  (ϵ/2)|(x/a)((|)^k)  the exact scaling relation applies:  (Subscript[ℰ, nk]) ∝ ϵ (η^( μ)) with η =(ℏ/(a(Sqrt[m ϵ]))) and μ=((2k)/(2+k))
If you know eigenvalue (Subscript[ℰ, nk]) for one set of parameters {ℏ, m, ϵ, a}, you can scale it exactly for any other values of {ℏ, m, ϵ, a}.
These estimates include the effect of the inverse mass parameter in the Schrödinger equation above.

If[homogeneous&& maybek>0 && !onesidedpotential,(
cc=.69;κ1=135;κ2=1427;(*from book appendix A - estimated ν parameter*);
νest=1/2(cc Exp[(2-k)/κ1]+(1-cc) Exp[(2-k)/κ2]) (*better for moderate n*);ν=νest;semiclassical=Table[ℰgb[n]/.psoln/.{ℏ->1,m->1/γ,ϵ->1}//N,{n,bslist}];(*the ϵ->1 is for cases in which the expression doesn't completely simplify*)Row[{Framed@ListPlot[{eigs,semiclassical},ImageSize->7*72,PlotMarkers->{"OpenMarkers",boldface@"+"},PlotLabel->"Phase Method values (O) vs Appendix A.3 semiclassical estimates(+) for ν ->"<>ToString[ν],Frame->True,GridLines->Full,FrameLabel->{"n","(Subscript[ℰ, n])"}],Framed@nf[ngooddigits]@TableForm[Transpose[{eigs,semiclassical}],TableHeadings->{None,{"PM","ℰgb(ν->"<>ToString[TraditionalForm@ν]<>")"}}]}])]

Solve for the wavefunctions if desired and do various checks

If[skipwavefunctions,(totalexecutiontime=(AbsoluteTime[]-totalstarttime);Print[Framed@("totalexecutiontime="<>ToString[totalexecutiontime/60 min])]; Print[Framed@("time/eigenvalue="<>ToString[totalexecutiontime/neigs sec])];allokbeep)]

If[skipwavefunctions,(SetOptions[#,CellOpen->True]&/@inputcells;SetOptions[#,CellOpen->True]&/@textcells;
Print[Framed@Framed@boldface@"skipwavefunctions flag is set to True, quitting kernel now"];Quit[ ])];

Print["{U, emaxsearch, tpsf, turningpoints[emaxsearch], {xmin,xmax}}=\n",{U,emaxsearch,tpsf,turningpoints[emaxsearch],{xmin,xmax}}//N];

this generates an estimate of the wavefunction which is normalized later. The initial slope/value doesn’t matter much at all, for reasons explained in the book.

Clear[getwf]

Print["ODE solver=",integrationmethod]

getwf[ℰ_]:=getwf[ℰ]=(zeros={};
soln=NDSolve[Schreq[ℰ],Φ[x],{x,xmin,xmax},Method->integrationmethod,WorkingPrecision->workingprecision,MaxSteps->10^7][[1]]//Quiet;(Φ[x]/.soln));

We will estimate the point at which the wavefunction starts to diverge at large x by finding the minimum of Log@Abs@wf, using a starting search point near xmax.
Here is an example.  The cutoff location in x is most sensitive to the accuracy of the eigenvalue, and hardly at all by the initial slope/value ratio. 
The cutoff slowly moves to higher x for higher precision eigenvalues.

nuthin=10^(-10);(*note: below, by displacing the lower turning point graphic upward by tiny amount 'nuthin' we avoid graphics error for some one-sided bottomless potentials at Machine Precision*)

Clear[xcutoff,cutoffs,pltcut,uwf,tmpfun];xcutoff[n_]:=xcutoff[n]=x/.(FindMinimum[{Evaluate[Log@Abs@getwf[eigs[[n]]]],x>Max[turningpoints[eigs[[n]]]]},{x,tplist[[2]],xmax*.99},Method->"InteriorPoint",AccuracyGoal->10]//Quiet)[[2]];
(*xcutoff[n_]:=xcutoff[n]=(tmpfun[x_]:=Log@Abs@uwf[n,x];xxmin=turningpoints[eigs[[n]]][[-1]];xx=xmax//N;δxx=.5;
TimeConstrained[While[(tmpfun[xx-δxx]<tmpfun[xx])&& (xx>=xxmin),xx=xx-δxx],1];Return[Max[xx,.1]])//Quiet*)
pltcut[n_]:=(xcut=xcutoff[n];uwf[x_]=getwf[eigs[[n]]];tplist=MinMax[turningpoints[eigs[[n]]]];
Plot[Chop@Log@Abs@Evaluate@uwf[x],{x,Max[xmin,tplist[[1]]-.5(tplist[[2]]-tplist[[1]])],Min[xmax,tplist[[2]]+3*(xcut-tplist[[2]])]},PlotRange->All,PlotPoints->1000,Frame->True,FrameLabel->{{"log(|Φ(x)|)",None},{"x","Φ(x) large x cutoff, n="<>ToString[n]<>", ℰ="<>ToString[NumberForm[eigs[[n]],{10,5}]]}},Epilog->{PointSize[.02],Point[{{tplist[[1]],Log[Abs[uwf[x]/.x->nuthin+tplist[[1]]]]},{tplist[[2]],Log[Abs[uwf[x]/.x->tplist[[2]]]]}}],InfiniteLine[{{xcut,1},{xcut,-1}}]}]);

This shows the estimated cutoffs for each of the eigenfunctions.  They are plotted so the user can verify that nothing went wrong. Sometimes it will.

ncols=Nearest[Divisors[neigs],4]//Last ;(*ncols is number of columns in which to display wavefunction plots; here calculate reasonable one*);
If[Not[3<=ncols<=5],ncols=4];(*in case of awkward column numbers, just pad it*)
placeholder=Graphics[{Gray,PointSize[.1],Point[{0,0}]}];(*graphic placeholder for padding elements in grid*)
allwfcutplot=GraphicsGrid[Partition[PadRight[ParallelTable[pltcut[n],{n,1,Length[eigs]}],ncols*Ceiling[neigs/ncols],placeholder],ncols],ImageSize->14*72]

This plot compares the cutoff locations with the upper turning point location. 
The cutoffs were set at the point at which the exponential runaway behavior starts, above the upper turning point. 
We see the cutoffs pretty much track with the upper turning points.

Clear[cutoffs];
cutoffs=ParallelTable[xcutoff[n],{n,1,Length[eigs]}];
uppertps=ParallelTable[turningpoints[eigs[[n]]][[-1]],{n,1,Length[eigs]}];
GraphicsRow[{ListPlot[{cutoffs,uppertps},Frame->True,PlotRange->All],ListLogLogPlot[{cutoffs/uppertps},Frame->True,PlotRange->All]//Quiet}]

Normalize the wavefunction to unity over the range the wavefunction solution is valid.   Outside that region, where it starts to diverge, define the wavefunction (as a Piecewise function) to be identically zero, so we don’t have to worry about the valid ranges from then on.  If desired this could be improved slightly by splicing onto a extrapolated exponential decay, but that is not done here, the values are small enough. Here the wavefunctions are “memo-ized” (cached) so they will be generated once, then can be used as needed.

AbsoluteTiming[Print["wf cutoffs=\n",ParallelTable[xcutoff[n],{n,1,neigs}]];]

Clear[wavefunction,norm];
wavefunction[n_]:=wavefunction[n]=(
xcut=xcutoff[n];
Piecewise[
{{0., x<xmin},{0., x>Min[xmax,xcut]}},
(uwf=Evaluate@getwf[eigs[[n]]];
norm=NIntegrate[uwf^2,{x,xmin,xcut},MaxRecursion->100,WorkingPrecision->workingprecision];
uwf/Sqrt[norm])])

Tabulate the normalized wavefunctions (as interpolated function objects) so they can used later
They are also saved as an Iconized object so they can be copied into another notebook and used there. 
You can use a syntax like this to evaluate on whatever x-grid you like.  	
	waves[[nwf,1]] gives the allowed range
	waves=theiconized object;nwf=1;
ParallelTable[{x,Evaluate[waves[[nwf]]]},{x,.1,10,.1}]

neigs=Length[eigs];AbsoluteTiming[allwfs=Table[wavefunction[n]//Quiet,{n,1,neigs}];]

Iconize[bslist](*these are the corresponding state numbers*)

Iconize[allwfs](*these can be copied to another notebook and evaluated on whatever suitable grid you like*)

Plot the lowest and highest energy wavefunctions
change the plot scale to whatever you want

{xplotmin,xplotmax}=3*tplist;

GraphicsRow[{Plot[wavefunction[1],{x,xplotmin,xplotmax},PlotRange->All],Plot[wavefunction[neigs],{x,xplotmin,xplotmax},PlotRange->All]}]

Check wavefunction orthonormality
You can turn off this by setting the flag noortho=True
Significant non-orthonormality is a sign of a problem; consider changing workingprecision if cutoff plots show problems. 
All it takes is for one wavefunction to be off. 
Conversely, if wavefunctions are orthonormal to a high degree it’s a sign the calculation has gone well.

Print["noortho=",noortho]

nintegrals=neigs*(neigs+1)/2;

Print["# of orthogonality integrals to do is:  ",nintegrals]

maxorthointegrals=500;(*intercept situations that would likely take an intolerably long time to run, in that case turn off orthogonality integrals; max value is up to user*)

If[nintegrals>maxorthointegrals,Print["Resetting noortho flag: too many ortho integrals to do in reasonable time:"]; noortho=True]

If[!noortho,(Print["orthonormality table timing (sec)=",orthotime=AbsoluteTiming[orthotable=Table[NIntegrate[wavefunction[i]wavefunction[j],{x,xmin,xmax},MaxRecursion->Automatic,WorkingPrecision->workingprecision]//Quiet,{j,1,neigs},{i,1,j}];]];
normerrmat=(Normal[LowerTriangularMatrix[orthotable]]-IdentityMatrix[neigs]);
Print["Average orthonormality error=",normerr=Norm[normerrmat]/(neigs*(neigs-1)/2),"  \nMinMax error=",ScientificForm[MinMax[Flatten[normerrmat]],{4,3}]];
Print[Framed@boldface@"Wave Function Orthonormality Table"];
Print[TableForm@N@Round[orthotable,10^(-3)]];
If[normerr>10^(-3),
Print[Framed@Framed@boldface@"Warning: the orthonormality matrix should be very close to unit matrix but it's not.  Check cutoff plots. Some wave functions are invalid or there may be duplicates."];
];),Print["orthonormality table calculation disabled by noortho flag"],Print["warning: noortho flag is not set to True or False"]]

here are normalized wavefunctions

Print["{U,neigs,xplotmin,xplotmax}=",{U,neigs,xplotmin,xplotmax}]

If the first eigenfunction found is supposed to be the ground state, and the potential does not have a barrier at the origin,
make sure the eigenfunction has no internal zeros (nodes), otherwise it’s not the ground state

Print["{bottomless, onesidedpotential}=",{bottomless,onesidedpotential}]

xbar=NIntegrate[x*wavefunction[1]^2,{x,xmin,xmax},AccuracyGoal->5];

If[bottomless,Plot[wavefunction[1],{x,xplotmin,3*xbar},PlotRange->All]]

change the plot range to whatever you want to

If[!(NumericQ[xplotmin]&&NumericQ[xplotmax]),{xplotmin,xplotmax}=3*tplist]

ncols=Nearest[Divisors[neigs],4]//Last ;(*ncols is number of columns in which to display wavefunction plots; here calculate reasonable one*);
If[Not[3<=ncols<=5],ncols=4];(*in case of awkward column numbers, just pad it*)
placeholder=Graphics[{Gray,PointSize[.1],Point[{0,0}]}];(*graphic placeholder for padding elements in grid*)
pltlbl[n_]:="n="<>ToString[n]<>",  ℰ="<>ToString[NumberForm[eigs[[n]],6]]<>",\n   xcutoff="<>ToString[NumberForm[xcutoff[n],5]]; 
wavefunctionplot=GraphicsGrid[Partition[PadRight[Table[(tplist=MinMax[turningpoints[eigs[[n]]]];Plot[wavefunction[n],{x,xplotmin,xplotmax},FrameLabel->{{"ψ(x)",None},{"x",pltlbl[n]}},PlotRange->All,Frame->True,PlotPoints->200,Epilog->{PointSize[.015],Point[{{tplist[[1]],0},{tplist[[2]],0}}]}]),{n,1,Length[eigs]}],ncols*Ceiling[neigs/ncols],placeholder],ncols],ImageSize->12*72]

this is an internal consistency check: are the expectation values of the Hamiltonian in each state in agreement with the computed eigenvalues?

Print["(ℏ^2)/m=\"γ\"=",γ];

Clear[H,ψ,energyexpectationvalue];(*values calculated below are cached for later use*)
ψ=wavefunction;Clear[energyexpectationvalue];
H:=(-γ/2 D[Evaluate[#],{x,2}]+U #)&; 
energyexpectationvalue[n_]:=energyexpectationvalue[n]=NIntegrate[ψ[n]H[ψ[n]],{x,xmin,xmax},MinRecursion->Automatic,MaxRecursion->50,WorkingPrecision->workingprecision,AccuracyGoal->ngooddigits+2]//Quiet

The following table is a comparison of the computed eigenvalues to the expectation values of the Hamiltonian which are calculated using the computed wave functions.  The discrepancy between them should not be taken as the accuracy of the eigenvalue. 
The computed eigenvalues in general will be considerably more accurate than the computed expectation values, which are just computed and shown here as a check in consistency between computed wavefunctions and corresponding eigenvalues

timelimit=240;
TimeConstrained[({timing,evs}=AbsoluteTiming[ParallelMap[energyexpectationvalue,Range[1,neigs]]];Print["timing=",timing,sec]; nf[ngooddigits]@TableForm[Transpose[{evs,eigs,evs/eigs-1 }],TableHeadings->{None,{"     <ℰ>","computed evals","relative discrepancy"}}]),timelimit,"Abort: this is taking too long"]

for homogeneous potential, checks that virial theorem and  <T>+<U> = ℰ are satisfied:
i.e. 2 <T> /( k <U>)=1 and (<T> + <U>)/ℰ = 1

Clear[ψ,KEbar,PEbar,virialratio];(*values calculated below are cached for later use*)

AbsoluteTiming[
If[homogeneous,
(k=maybek;ψ=wavefunction;Toper:=(-γ/2 D[Evaluate[#],{x,2}])&; Uoper:=(U #)&; 
KEbar[n_]:=KEbar[n]=NIntegrate[ψ[n]Toper[ψ[n]],{x,xmin,xmax},MinRecursion->Automatic,MaxRecursion->Automatic,WorkingPrecision->workingprecision,AccuracyGoal->Max[8,ngooddigits]]//Quiet;
PEbar[n_]:=PEbar[n]=NIntegrate[ψ[n]Uoper[ψ[n]],{x,xmin,xmax},MinRecursion->Automatic,MaxRecursion->Automatic,WorkingPrecision->workingprecision,AccuracyGoal->Max[8,ngooddigits]]//Quiet;
virialratio[n_]:=virialratio[n]=(2*KEbar[n])/(k*PEbar[n]);
Print[separator];
Print[Framed[boldface@"the following values, evaluated for each wavefunction, should all be close to 1.0"]];
Print[boldface@"virial ratios=\n",nf[6]@Table[virialratio[n],{n,1,neigs}]];Print[boldface@"energy total ratios (<T>+<U>)/ℰ=\n",nf[6]@(Table[KEbar[n]+PEbar[n],{n,1,neigs}]/eigs)];
)]]

if these relative errors are small, eigenvalues and expectation values of Hamiltonian in each wavefunctions are consistent

Next,  check whether  the range is sufficient to sufficiently converge the highest energy WF to zero; here we calculate the probability of the highest wavefunction if we shrink the range slightly.
I just test at the high end because, for one-sided potentials,  it’s inappropriate to check at xmin because WF is forced to go to zero there but may be substantial near x=0.

{xmin,xmax}

{xplotmin,xplotmax}

xbar=(xmin+xmax)/2;
highestwfnormerr=Norm[1-NIntegrate[wavefunction[neigs]^2,{x,xmin,xbar+0.95(xmax-xbar)},MinRecursion->Automatic,MaxRecursion->100,WorkingPrecision->workingprecision,AccuracyGoal->Automatic]//Quiet]

probability density for highest computed energy level

Plot[wavefunction[neigs]^2,{x,xmin,xmax},PlotRange->{{xplotmin,xplotmax},All},PlotPoints->100]

If[highestwfnormerr>10^(-7),warningbeep;
Print[Framed@Framed@"Warning: x-range may not be sufficient to ensure highest WF convergence - try extending the range (increase upper tpsf)"]]

Finally, as a measure of quantum tunneling (to correlate with the (Subscript[ν, kn]) parameter in book) calculate the total probability for each states that resides outside the classical turning points.  Here we evaluate this only for homogeneous potentials and at $MachinePrecision which should not be very time-consuming. It’s cached for reuse.

Clear[proboutside];

proboutside[n_]:=proboutside[n]=({x1,x2}=MinMax[turningpoints[eigs[[n]]]];1-NIntegrate[ψ[n]^2,{x,x1,x2},MinRecursion->Automatic,MaxRecursion->Automatic,WorkingPrecision->$MachinePrecision,AccuracyGoal->Max[8,ngooddigits]]//Quiet);

If[homogeneous,probsoutside=ParallelTable[proboutside[n],{n,1,neigs}];Print@ListPlot[probsoutside,Frame->True,GridLines->Full,PlotRange->All,FrameLabel->{"n","prob outside"},PlotLabel->"Probability vs n to find particle outside classical turning points"];
Print@TableForm[Transpose[{bslist,probsoutside}]]];

If[homogeneous,Print["and here are the iconized values of probability outside classical turning points"];Iconize[probsoutside]]

write out plots and data

filenametag=(ToString[U,InputForm]<>" "<>integrationmethod<>" WP"<>ToString[Round[workingprecision]])//StringReplace["/"->"_over_"](*this can be used to generate output file names*);

If[onesidedpotential,filenametag=filenametag<>" one-sided"]

Print["filenametag=",filenametag];

Export["__most recent potential.txt",filenametag];(*this just creates a text file containing the name of the potential*)

outputdirectory=filenametag<>" outputfiles";

If[writefiles,Print["outputdirectory=",outputdirectory]];

If[!FileExistsQ[outputdirectory]&& writefiles,CreateDirectory[outputdirectory]]

plottypes={potentialplot,phaseplt,siftplot,eigplot,levelplot,allwfcutplot,wavefunctionplot};

eigtbl=TeXForm@TableForm[Partition[nf[ngooddigits]/@eigs,10],TableDirections->Row];

If[writefiles,
(SetDirectory[outputdirectory];
schreqoutfile=Export[filenametag<>" schreq.pdf",TraditionalForm@{Schreq[ℰ],"\n",H[Φ[x]]},"PDF"];
potentialoutfile=Export[filenametag<>" potential.pdf",potentialplot,"PDF"];
phaseplotoutfile=Export[filenametag<>" phaseplot.pdf",phaseplt,"PDF"];
siftoutfile=Export[filenametag<>" sift.pdf",siftplot,"PDF"];
eigoutfile=Export[filenametag<>" eig.pdf",eigplot,"PDF"];
leveloutfile=Export[filenametag<>" level.pdf",levelplot,"PDF"];
allcutoutfile=Export[filenametag<>" allcut.pdf",allwfcutplot,"PDF"];
wavefunctionoutfile=Export[filenametag<>" wavefunction.pdf",wavefunctionplot,"PDF"];
eigslistoutfile=Export[filenametag<>" eiglist.txt",NumberForm[eigs,ngooddigits],"Table"];
eigstableoutfile=Export[filenametag<>" eigtableTeX.txt",eigtbl,"Lines"];
)]

these are the output file names for plots etc

FileNames[filenametag<>"*"]//ColumnForm

phaseplt (*"phaseplot" variable is the function, this is the plot of it*)

siftplot

eigplot

write output files (if writefiles flag is True)

If[writefiles,(Export[filenametag<>" phaseplot.eps",phaseplt,"EPS"];
Export[filenametag<>" sift.eps",siftplot,"EPS"];
Export[filenametag<>" eig.eps",eigplot,"EPS"];)]

Print["{U, neigs, eminsearch, emaxsearch, xmin, xmax}=\n",{U,neigs,eminsearch,emaxsearch,xmin,xmax}]

Finish and clean up

Print time stamp, total execution time, do cleanup, save eigenvalues as an Iconized object

{sifttime,refinetime,sifttime+refinetime,orthotime[[1]]}/60 minutes

totalexecutiontime=(AbsoluteTime[]-totalstarttime);

Framed[totalexecutiontime/3600 hours]

Framed[totalexecutiontime/60 minutes]

Framed[(totalexecutiontime/neigs) seconds/eigenvalue]

Iconize[eigs]

Print["end-of-execution date/time: ",Now]

elapsedtime=Round[(AbsoluteTime[]-totalstarttime)/60.,.1] minutes;

Print[Framed@("total elapsed time= "<>ToString[elapsedtime])];
{$Version,$OperatingSystem,$ProcessorType,$ProcessorCount};
speakthis="total elapsed time was "<>ToString[elapsedtime];
allokbeep;
If[speakelapsedtime,AudioPlay[SpeechSynthesize[speakthis]];]
(*If[speakelapsedtime,Speak[speakthis]]*)

Do some cleanup, close palette, or not

(*restore all code and text cells to their open states, so they're accessible once notebook is closed and then reopened, and delete palette so they don't pile up*)

Print[Framed@Style["Note: if text and code cells were left closed, executing the following (very last) cell below\!\(\*StyleBox[\(\\ \*StyleBox[" ",FontWeight->"Bold"]\)]\)⇓ will open them",Bold]];

nb=SelectedNotebook[];
inputcells=Cells[nb,CellStyle->"Input"];
textcells=Cells[nb,CellStyle->"Text"];
SetOptions[#,CellOpen->True]&/@inputcells;SetOptions[#,CellOpen->True]&/@textcells;

closeit=ChoiceDialog["Close Navigation Palette?  \[Implies]",{"Close It"->True,"Keep it Open"->False},WindowFloating->True,WindowMargins->{{.65*screenwidth,Automatic},{Automatic,0}}];
If[closeit,Do[NotebookClose[navigationpalette],{3}]];