Setting some parameters

ClearAll[R,Mn,Edeut,Vzero]
ϵ=N[(10^(-4))];
hbarc=198;
R=1.93/hbarc;
Mn=940.;
Edeut=2.225;

Solving transcendental equation to get the square well height in MeV

Vzero=Vzero/.FindRoot[Cot[(Sqrt[Mn(Vzero-Edeut)])R]==-(Sqrt[(Edeut/(Vzero-Edeut))]),{Vzero,20}]

Defining the square well potential

Vtotal[r_]=If[r<R,-Vzero,0];
potplot=Plot[Vtotal[r],{r,0,3R},PlotStyle-> Red]

Solving the Schrodinger equation numerically for the bound state

normali=0.75 (10^4);     (* Normalization factor for the wavefunction *)
crap=NDSolve[{(-(1/Mn)u''[r]+(Vtotal[r]+Edeut)u[r])==0,u[ϵ]==ϵ,u'[ϵ]==1},u,{r,ϵ,5R}];
BSplot=Plot[Evaluate[normali(u[r]/.crap)],{r,ϵ,5R},Frame->True,PlotStyle->Green]

Scattering states: Finding the scattering length

Let' s first plot the wavefunction at zero energy

k=0;
normali=0.75 (10^4);    (* Normalization factor for the wavefunction *)
crap=NDSolve[{(-(1/Mn)u''[r]+(Vtotal[r]-(k^2))u[r])==0,u[ϵ]==ϵ,u'[ϵ]==1},u,{r,ϵ,20R}];
scattplot=Plot[Evaluate[normali(u[r]/.crap)],{r,ϵ,5R},Frame->True,PlotStyle->Blue]

Now let' s get the scattering length using method described in class

ClearAll[crapA,rtest,Deriv,atest,ascatt];
ascatt:=Module[{atest},
crapA=NDSolve[{(-(1/Mn)u''[r]+Vtotal[r]u[r])==0,u[ϵ]==ϵ,u'[ϵ]==1},u,{r,ϵ,20R}];
rtest=4R;
Deriv= D[(u[r]/.crapA),r]/.r->rtest;
atest=(rtest-((u[r]/.crapA)/.r->rtest)/Deriv)]

Recall that experimentally a = 5.4 fm

ascatt[[1]] hbarc

Combined plot

Show[scattplot,BSplot,potplot]

The nucleons in the deuteron live most of the time outside of the range of the potential that binds them!