Find the bound state energy levels (i.e., stationary states
with E<0) of a proton in a finite square well with depth=50MeV,
full width=2a=10fm.   Begin by finding the dimensionless quantity z_0,
using  Mathematica to solve Eqn 2.156 and the equivalent equation
for the odd states, and then converting the resulting z
to energy (in MeV).  Having Mathematica plot both sides
of the equation to be solved (with FindRoot[]), allows you to
to see the solutions and thus give Mathematica a better initial
guess at a solution.

z0 is essentially a dimensionless measure of the depth of potential

z0= 5e-15*sqrt(2*mp*50e6*e)/hbar = 7.7615

Note: used online calculator from www.physics.csbsju.edu;
fm="fermi"=1e-15 m 
e is the conversion constant between eV and J

z0=7.7615
Plot[{Tan[z],Sqrt[(z0/z)^2-1]},{z,.2,10}]

see intersection at z approx 1.4, 4, 6.7
FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,1.4}]
Out[3]= {z -> 1.39065}

FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,4}]
Out[4]= {z -> 4.14848}

FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,6.7}]
Out[5]= {z -> 6.78913}


for odd states:
Plot[{Cot[z],-Sqrt[(z0/z)^2-1]},{z,.2,10}]
see intersections: 2.7, 5.5

FindRoot[Cot[z]==-Sqrt[(z0/z)^2-1],{z,2.7}]
Out[11]= {z -> 2.77585}

FindRoot[Cot[z]==-Sqrt[(z0/z)^2-1],{z,5.5}]

Out[12]= {z -> 5.49631}

allZ={Last[Last[%3]],Last[Last[%11]],Last[Last[%4]],Last[Last[%12]],Last[Last[%5]]}

Out[15]= {1.39065, 2.77585, 4.14848, 5.49631, 6.78913}


Energy = (hbar/a)/(2 m) z^2 - 50MeV
seek: (hbar/a)/(2 m) in usits of MeV:

(hbar/5e-15)^2/(2*mp)/(e*1e6)
0.82999 

%15^2*0.82999-{50,50,50,50,50}
Out[16]= {-48.3949, -43.6046, -35.716, -24.9265, -11.7439}

