This notebook implements what might be called the "Shoot First" method, a kind of shooting method, with bisection, to home in on a single eigenvalue.
Like the Phase Method (here we borrow some approaches from the PM such as turning point scale factors).
 It works by starting the initial value problem so far into the classically forbidden region
that you don’t have to worry about the initial conditions (for reason explained in the Phase Method Book);
 just let it evolve, and evaluate the sign of the solution at xmax, and determine whether sign flips occur between 
 endpoints of a specified energy interval that is intended to be smaller than the
spacing between the levels.  Bisection is then used to home in a specific eigenvalue.
Unlike the PM this does not use the Phase Plot to count the states below a given energy, it just uses the sign flip
(Griffiths’s “wag the dog”). But in contrast with most other shooting method approaches (including Griffiths’s “wag the dog”) , 
no effort is made to get the solution to satisfy the boundary conditions or to guess what they might be. 
Whence the “Shoot First” (“Shoot first and ask questions later”). 

This works fine if the specified energy interval is known to contain only a single eigenvalue.
If not, an even number of eigenvalues in an interval will look like zero eigenvalues, and an odd number of eigenvalues looks like a single one.

In such cases it’s better to  use Phase Method which maps them all out automatically without ambiguity.

For well-behaved potentials ImplicitRungeKutta gives very good results at machine precision but it’s brittle and can give bad answers 
for other potentials e.g. Coulomb. LSODA at extended precision is quite robust.  Consider yourself warned.

Here change the potential to whatever you want 
but make sure the range extends well past turning points for highest energy chosen
(this is achieved here using turningpoint scale factors)

starttime=AbsoluteTime[ ];

U=x^2/2;
workingprecision=25 ;(*$MachinePrecision is fastest but accuracy is lower*)
SP=SetPrecision[#,workingprecision]&;
integrationmethod="LSODA";(*LSODA, BDF, or if potential not singular, ImplicitRungeKutta*)

etest=10;(*upper energy limit to find eigenvalues*)

turningpoints=MinMax[x/.NSolve[U==etest,x,Reals]]

tpsf={3,3};(*turning point scale factors*)

{xmin,xmax}=tpsf*turningpoints

Plot[{U,etest},{x,xmin,xmax}]

The following code memo-izes (caches) previously computed values of Φ(xmax) to avoid redundant computations. 
This speeds up the bisection process nearly four-fold.

Clear[Phiofxmax];
ℰ=.;Φ0p=.;Φ0=1/10;
Phiofxmax[ℰ_,Φ0p_]:=Phiofxmax[ℰ,Φ0p]=(soln=Quiet@NDSolve[{Φ''[x]+2 (ℰ-U)Φ[x]==0,Φ[xmin]==Φ0,Φ'[xmin]==Φ0p},Φ[x],{x,xmin,xmax},Method->integrationmethod,WorkingPrecision->workingprecision,MaxSteps->10^7][[1]];Return[Φ[x]/.soln/.x->xmax])

signflip[{ℰ1_,ℰ2_}]:=(Sign[Phiofxmax[ℰ1,1/10]]*Sign[Phiofxmax[ℰ2,1/10]]==-1);
bisect[{ℰ1_,ℰ2_}]:=(ebar=(ℰ1+ℰ2)/2;If[signflip[{ℰ1,ℰ2}],If[signflip[{ℰ1,ebar}],{ℰ1,ebar},{ebar,ℰ2}],Throw["zero, or an even number of eigenvalues are in interval\nor the x-range is incorrect"]])

eps=SP[.0001];(*for extending interval a hair to catch eigenvalue that falls right on interval boundary*)
(*nits=3+Log[2,10]*workingprecision/2//Ceiling;(*with lsoda ultimate accuracy seems to be about workingprecision/2*)*)
(*geteig[{e1_,e2_}]:=(Catch[Mean@Nest[bisect,{e1-eps,e2+eps},nits]])*)
geteig[{e1_,e2_}]:=(Catch[Mean@FixedPoint[bisect,{e1-eps,e2+eps}]])

here is a simple example - find the single eigenvalue presumed to lie between 8 and 9

answer=AbsoluteTiming[geteig[{8,9}]]

integrationmethod

If[integrationmethod=="LSODA",ngooddigits=Floor[(workingprecision-1)/2],ngooddigits=Round@workingprecision];
If [integrationmethod!="LSODA", Print[Framed["Caveat Emptor for ngooddigits when integrationmethod = "<>integrationmethod]]];

ngooddigits

here’s your one answer

NumberForm[answer[[2]],ngooddigits]

All of the following stuff is for illustrative purposes and can be deleted if you just want to use the executable to find eigenvalues.
The comparison with exact values is specific to harmonic oscillator U=x^2/2.
This is a brute-force search for all the eigenvalues of an harmonic oscillator, from 0 to etest,  in many small-ish search intervals.
This is inefficient but suffice if the eigenvalues are spaced out well enough and the search intervals are fine enough.
Otherwise the Phase Method is much better for resolving them all automatically, efficiently, and systematically.
See PM notebooks for that.

fuzz=.001Random[];

searchintervals=SP@Table[{e,e+.1},{e,0,etest,.1-fuzz}];

Length[searchintervals]

AbsoluteTiming[results=ParallelMap[geteig,searchintervals];]

results//Length

here are the found eigenvalues

eigs=Select[results,NumericQ]

exact=Table[n-1/2,{n,1,Length[eigs]}];

relative error compared to exact

eigs/exact-1

ListLogPlot[Abs[eigs/exact-1],Frame->True,GridLines->Full,PlotRange->All,FrameLabel->{"n","relative error"}]

roundedeigs=N[eigs,ngooddigits]

Precision[roundedeigs]

Iconize[roundedeigs]

elapsedtime=(AbsoluteTime[]-starttime)/60. min

The following stuff shows that the initial conditions don’t matter as long as you start the ODE solution well below the first turning point.
The values of the solutions get into the stratosphere which is a good thing.  The mathematical reason can be found in the PM book. 
This section of this notebook is illustrative and can be deleted if you just want to find eigenvalues.

This demonstrates the very rapid convergence of the logarithmic derivatives for two different initial slopes
as long as the integration is started well out into the classically forbidden region (below the first turning point).
Initial conditions don’t really matter under such conditions.

Φfun[ℰ_,Φ0p_]:=(soln=Quiet@NDSolve[{Φ''[x]+2 (ℰ-U)Φ[x]==0,Φ[xmin]==1/10,Φ'[xmin]==Φ0p},Φ[x],{x,xmin,xmax},Method->integrationmethod,WorkingPrecision->workingprecision,MaxSteps->10^7][[1]];Return[Φ[x]/.soln])

Λ[ℰ_,Φ0p_]:=D[Log@Φfun[ℰ,Φ0p],x]

LogPlot[Evaluate@{Λ[5/2,10]/Λ[5/2,.1]},{x,xmin,xmin+1},PlotRange->All]

Equivalently, this shows the solutions Φ themselves very rapidly (within Δx ~ 0.1) approach a constant ratio under those conditions. 
The constant ratio is irrelevant -- it goes away when wavefunctions are normalized.

LogPlot[Evaluate@{Abs@Φfun[5/2,.1],Abs@Φfun[5/2,10]},{x,xmin,xmin+1}]

LogPlot[Evaluate@{Abs@Φfun[5/2,.1],Abs@Φfun[5/2,10]},{x,0,xmax}]

MinMax@ParallelTable[Phiofxmax[ℰ,1/10],{ℰ,fuzz,10,.1+fuzz}]