(* Mathematica program for Problem 8.2

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

R = Function[{n, l, Z, r}, (2 Z/(n a0))^(3/2) *
	     Sqrt[(n - l - 1)!/(2 n)/((n + l)!) ] (2 Z r/(n a0))^l *
            LaguerreL[n - l - 1, 2 l + 1, 2 Z r/(n a0)] E^(-Z r/(n a0))]

(* Radial Functions *)

Print["R10(r)=",R[1,0,1,r]]
Print["R20(r)=",R[2,0,1,r]]
Print["R21(r)=",R[2,1,1,r]]

(* Plots of Rnl (r) *)
Plot[R[1,0,1,r] /. a0->1,{r,0,10}, PlotRange->{0,2},
AxesLabel->{"r/a0","R10(r)"}]
Plot[R[2,0,1,r] /. a0->1,{r,0,10},PlotRange->{-1,1},
AxesLabel->{"r/a0","R20(r)"}]
Plot[R[2,1,1,r] /. a0->1,{r,0,10},PlotRange->{0,0.2},
AxesLabel->{"r/a0","R21(r)"}]

(* Plots of Radial Probability *)
Plot[r^2 R[1,0,1,r]^2 /. a0->1,{r,0,10}, PlotRange->{0,0.6}, 
AxesLabel->{"r/a0","r^2 R10(r)"}]
Plot[r^2 R[2,0,1,r]^2 /. a0->1,{r,0,10},PlotRange->{0,0.2},
AxesLabel->{"r/a0","r^2 R20(r)"}]
Plot[r^2 R[2,1,1,r]^2 /. a0->1,{r,0,10},PlotRange->{0,0.2},
AxesLabel->{"r/a0","r^2 R21(r)"}]
Plot[r^2 R[3,0,1,r]^2 /. a0->1,{r,0,30},PlotRange->{0,0.12},
AxesLabel->{"r/a0","r^2 R30(r)"}]

(* Normalization *)
Print["Integral r^2 R10(r)^2 = ",
      Integrate[r^2 R[1,0,1,r]^2,{r,0,Infinity},Assumptions->a0>0]]
Print["Integral r^2 R20(r)^2 = ",
      Integrate[r^2 R[2,0,1,r]^2,{r,0,Infinity},Assumptions->a0>0]]
Print["Integral r^2 R21(r)^2 = ",
      Integrate[r^2 R[2,1,1,r]^2,{r,0,Infinity},Assumptions->a0>0]]

(* Orthogonality *)
Print["Integral r^2 R10(r) R20(r) = ",
      Integrate[r^2 R[1,0,1,r] R[2,0,1,r],{r,0,Infinity},Assumptions->a0>0]]
Print["Integral r^2 R21(r) R31(r) = ",
      Integrate[r^2 R[2,1,1,r] R[3,1,1,r],{r,0,Infinity},Assumptions->a0>0]]
Print["Integral r^2 R10(r) R21(r) = ",
      Integrate[r^2 R[1,0,1,r] R[2,1,1,r],{r,0,Infinity},Assumptions->a0>0]]