(* Mathematica program for Problem 11.3

  date : 03/12/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(* Definition od Constants *)
alfa=1./137;
cnum=3*10^8;
munum=0.51*10^6/cnum^2;
a0num=197.35*10^(-9)/(munum*cnum^2*alfa);
Rnum=10.^(-15);
H0num=alfa munum cnum^2;

(* Definiton of the Radial Functions *)
Rnl = Function[{n, l, r}, (2/(n a0))^(3/2) Sqrt[(n - l - 1)!/(2 n)/((n + l)!) ]
 (2 r/(n a0))^l  LaguerreL[n - l - 1, 2 l + 1, 2 r/(n a0)] E^(-r/(n a0))];

(* Perturbation *)
Hp=H0 (-3/2 a0/R+ a0/r + 1/2 a0*r^2/R^3);

ans10=Integrate[r^2 Rnl[1,0,r]^2 Hp,{r,0,R}];
ans20=Integrate[r^2 Rnl[2,0,r]^2 Hp,{r,0,R}];
ans21=Integrate[r^2 Rnl[2,1,r]^2 Hp,{r,0,R}];

aux10= ans10 /. R->y a0;
aux10= Normal[Series[aux10,{y,0,2}]];
aux10= aux10 /y^2  /. {H0->μ c^2 α} ;
aux20= ans20 /. R->y a0;
aux20= Normal[Series[aux20,{y,0,4}]];
aux20= aux20 /y^4  /. {H0->μ c^2 α} ;
aux21= ans21 /. R->y a0;
aux21= Normal[Series[aux21,{y,0,4}]];
aux21= aux21 /y^4  /. {H0->μ c^2 α} ;

Print["Δ E10 = ",aux10," (R/a0)^2"]
Print["Δ E20 = ",aux20," (R/a0)^4"]
Print["Δ E21 = ",aux21," (R/a0)^4"]

num=SetPrecision[{a0->a0num,R->Rnum,H0->H0num,c->cnum},45];

res10= ans10 /. num ;
res20= ans20 /. num ;
res21= ans21 /. num ;

Print["Numerical Values: "]
Print["Δ E10 =",res10]
Print["Δ E20 =",res20]
Print["Δ E21 =",res21]