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=a*sqrt(2*mp*V*c^2)/(hbar*c)

on-line constant calculator:
Conversion constant: hbar*c/e = 197.33 MeV·fm 
5*sqrt(2*mp2*50)/197.33
z0=7.7614 

or using MKS units:
5e-15*sqrt(2*mp*50e6*e)/hbar
7.7615

even:
Plot[{Tan[z],Sqrt[(z0/z)^2-1]},{z,0,10}]
see intersections near: 6.7, 4.2, 1.4

FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,6.7}]
Out[43]= {z -> 6.78911}
FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,4.2}]
Out[44]= {z -> 4.14847}
FindRoot[Tan[z]==Sqrt[(z0/z)^2-1],{z,1.4}]
Out[45]= {z -> 1.39065}

odd:
Plot[{-Cot[z],Sqrt[(z0/z)^2-1]},{z,0,10}]
see intersections: 5.5,2.7
note that there is a missed intersection near z0;
if z0=2.5 Pi=7.85 we would have had one more state

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

now (ka)^2 = z0^2-z^2 and
E=-(hbar k)^2/(2m) = -hbar^2/(2*mp*a^2)*(z0^2-z^2)
= -(hbar*c)^2/(2*mp*c^2*a^2)*(z0^2-z^2)

online calculator:
-197.33^2/(2*mp2*5^2)*(7.7615^2-x^2)

f( 6.7891 )= -11.744
f( 5.4963 )= -24.927
f( 4.1485 )= -35.717
f( 2.7759 )= -43.605 
f( 1.3907 )= -48.396 
energies in MeV

alternatively use MKS:
-(hbar/5e-15)^2/(2*mp)*(7.7615^2-x^2)/e
f( 6.7891 )= -1.1743 × 10^7
f( 5.4963 )= -2.4926 × 10^7
f( 4.1485 )= -3.5715 × 10^7
f( 2.7759 )= -4.3604 × 10^7
f( 1.3907 )= -4.8394 × 10^7 
energies in eV

