u[r] was given in class12, R=u(r)/r
R[n_,l_,r_]=1/(n Sqrt[Pochhammer[(n-l), 2 l +1]]) (2 r/n)^(l+1) LaguerreL[n-l-1,2 l +1,2 r/n] Exp[-r/n]/r

Integrate[R[4,3,r] r^2 R[4,3,r] r^2, {r,0,Infinity}]Integrate[Integrate[Conjugate[SphericalHarmonicY[3, 3, theta, phi]] Cos[theta]^2 SphericalHarmonicY[3, 3, theta, phi] Sin[theta],{phi,0,2 Pi}],{theta,0,Pi}]
Out[2]= 40

Integrate[R[4,3,r]  R[4,3,r] r^2, {r,0,Infinity}]Integrate[Integrate[Conjugate[SphericalHarmonicY[3, 3, theta, phi]]  SphericalHarmonicY[3, 3, theta, phi] Sin[theta],{phi,0,2 Pi}],{theta,0,Pi}]

Out[3]= 1

Integrate[R[4,1,r] r^2 R[4,1,r] r^2, {r,0,Infinity}]Integrate[Integrate[Conjugate[SphericalHarmonicY[1, 0, theta, phi]] Cos[theta]^2 SphericalHarmonicY[1, 0, theta, phi] Sin[theta],{phi,0,2 Pi}],{theta,0,Pi}]
Out[4]= 360

length units are Bohr Radius (a), so for above e.g.,
<z^2> = 40 a^2
