
9.1:
from the past:

Recall the H-atom wavefunction in Mathematica form:
psi[n_,l_,m_,r_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]]) (2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n] SphericalHarmonicY[l, m, theta,phi]

it will help to have the complex conjugate of this wavefunction:
psic[n_,l_,m_,r_]=2/(n^2 Sqrt[Pochhammer[n-l,2l+1]]) (2r/n)^l LaguerreL[n-l-1,2l+1,2r/n] Exp[-r/n] SphericalHarmonicY[l, m, theta,-phi]

"r" here is in units of the Bohr radius (often denoted "a")
In general there are n^2 degenerate states for principal quantum number n

For n=2->4.  here is a list of the lm values:

lm={{0,0},{1,1},{1,0},{1,-1}}

QUESTION: using the sp encoding, the first one is 2s, the remaining three are 2p

9.1 asks you to find 4 integrals:
<2s|H|1s>, <2p,m=1|H|1s>, <2p,m=0|H|1s>, <2p,m=-1|H|1s>

where H = eEz...the eE is just an overall factor; if you want you can put that in later. z=r Cos[theta]
The number we get from the above wavefunction is in units of the Bohr radius


Table[
Integrate[Integrate[Integrate[psic[2,First[lm[[i]]],Last[lm[[i]]],r] r Cos[theta] psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
{i,1,4}]

               128 Sqrt[2]
Out[5]= {0, 0, -----------, 0}
                   243

if you don't like the fancy work to do these four integrals in one command, note that this is the same as:
Integrate[Integrate[Integrate[psic[2,0,0,r] r Cos[theta] psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
Integrate[Integrate[Integrate[psic[2,1,1,r] r Cos[theta] psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
Integrate[Integrate[Integrate[psic[2,1,0,r] r Cos[theta] psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
Integrate[Integrate[Integrate[psic[2,1,-1,r] r Cos[theta] psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],

9.11 starts the same: we need 4x3 integrals ({x,y,z} forming the vector r, and we need p= -e <j|vec r|i> the electric dipole)
<2s|x|1s>, <2s|y|1s>, <2s|z|1s>
<2p,m=1|x|1s>,<2p,m=1|y|1s>,<2p,m=1|z|1s>,
<2p,m=0|x|1s>,<2p,m=0|y|1s>,<2p,m=0|z|1s>,
<2p,m=-1|x|1s>,<2p,m=-1|y|1s>,<2p,m=-1|z|1s>,

Table[
Integrate[Integrate[Integrate[psic[2,First[lm[[i]]],Last[lm[[i]]],r] 
{ r Sin[theta] Cos[phi], r Sin[theta] Sin[phi] , r Cos[theta]} 
psi[1,0,0,r]  r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
{i,1,4}]

                       128   128 I             128 Sqrt[2]    128  128 I
Out[7]= {{0, 0, 0}, {-(---), -----, 0}, {0, 0, -----------}, {---, -----, 0}}
                       243    243                  243        243   243


The norm of all three 2p r vectors is 128 Sqrt[2]/243; the 2s r vector is zero...
Table[Norm[ %7[[i]] ], {i,1,4}]

             128 Sqrt[2]  128 Sqrt[2]  128 Sqrt[2]
Out[11]= {0, -----------, -----------, -----------}
                 243          243          243

the rest of the problem is basically unit conversion:

lifetime = 3 Pi epsilon0 hbar c^3/(omega^3 p^2)

where: omgea=(1/2) alpha^2 me c^2 (1-1/4)/hbar = (Delta E)/hbar
       p=e 128 Sqrt[2]/243 a

       plug into web site that know constants: http://www.physics.csbsju.edu/cgi-bin/twk/examples/F2.html
       
3*Pi*epsilon0*hbar*c^3/((.5*alpha^2*me*c^2*.75/hbar)^3*(e*128*Sqrt(2)/243*a )^2)
f( 0 )= 1.5953 × 10-9 sec