Michel Potential

Initialization

v0=21;
rr=5;
ρ=6.25;
rc=1.3;
z1=2;
z2=28;
e=1.6 (10^(-19));
αr=0.6;
α=18;
h=6.63 (10^(-34));
ℏ=h/π;
m=4;

Function Definition

vc[r1_]:=(((z1 z2 (e^2))/(2 rc)))(3-((r1^2)/(rc^2)))/;r1≤rc;
vc[r1_]:=((z1 z2 (e^2))/r1)/;r1>rc;
vr[r_]:=-v0 (1+α Exp[-((((r/ρ)))^2)])/(((1+Exp[((r-rr)/(2 αr))]))^2)+vc[r];

Ploting the Potential

Plot[vr[r],{r,0,10}]

rate=D[vr[r],r];
Plot[rate,{r,0,10}]

eα=14.47;
sln=NDSolve[{ψ''[r]+(eα-vr[r])ψ[r]==0,ψ[0]==0,ψ'[0]==1},ψ[r],{r,0,10}]
Plot[ψ[r]/.sln,{r,0,10},Axes->True]