(* Mathematica program for Problem 8.7

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

(* Definitions *)
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))];
f = Function[r, (a/Sqrt[Pi])^(3/2) E^(-a^2 r^2/2)];

(* Check Normalization *)
Print["Norm = ",Integrate[r^2 f[r]^2 4 Pi,{r,0,Infinity},Assumptions->a>0]]

(* Plot f and R10 for \alpha=a0=1 *)
Plot[{f[r] /. a->1,R[1,0,1,r]/Sqrt[4 Pi] /. a0->1},{r,0,5},AxesLabel->{"r/a0","f,R10"}]

(* Coefficient in ground state *)
A0=Integrate[f[r]*4*Pi*r^2*R[1,0,1,r]/Sqrt[4*Pi],{r,0,Infinity},
	     Assumptions->{a>0,a0>0}];
Print["A0 = ",A0]

(* Define \alpha = y/a0 and plot probablity as a function of y*)

A0Fun=Function[y,A0 /. a->y/a0];
Plot[A0Fun[y]^2,{y,0.1,10},PlotRange->{0,1},
     AxesLabel->{"y=α a0","|A0|^2"}]