we'll do 6.37 with Mathematica.
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]

In general there are n^2 degenerate states for principal quantum number n
For n=3->9.  here is a list of the lm values:

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

QUESTION: using the spd encoding, circle and label the 3s, three 3p, five 3d states.

Check the orthonormality of these states by calculating the 81 integrals: <n l' m'|n l m> where n=3
Do this in one chunk:
Table[Table[
Integrate[Integrate[Integrate[psic[3,First[lm[[i]]],Last[lm[[i]]],r]  psi[3,First[lm[[j]]],Last[lm[[j]]],r] r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
{i,1,9}],{j,1,9}]
MatrixForm[%]

QUESTION: explain (write down) why each term in the integral is the correct term.

For the Stark Effect we seek the z matrix elements:  <n l' m'|z|n l m>
Calculate these 81 integrals in one chunk:
Table[Table[
Integrate[Integrate[Integrate[psic[3,First[lm[[i]]],Last[lm[[i]]],r] r Cos[theta] psi[3,First[lm[[j]]],Last[lm[[j]]],r] r^2 Sin[theta],{phi,-Pi,Pi}],{theta,0,Pi}],{r,0,Infinity}],
{j,1,9}],{i,1,9}]

MatrixForm[%]
QUESTION: explain why each term in the integral is the correct term.

     0       0  -3 Sqrt[6]   0       0       0       0       0       0

     0       0       0       0       0     -9/2      0       0       0

-3 Sqrt[6]   0       0       0       0       0  -3 Sqrt[3]   0       0

     0       0       0       0       0       0       0     -9/2      0

     0       0       0       0       0       0       0       0       0

     0     -9/2      0       0       0       0       0       0       0

     0       0  -3 Sqrt[3]   0       0       0       0       0       0

     0       0       0     -9/2      0       0       0       0       0

     0       0       0       0       0       0       0       0       0

QUESTION: print out the above answer for the matrix.  Note that in the problem Griffiths gives
three check integrals.  Find and circle each one of these example integrals making sure
you get the order right...i.e., <A|z|B> is at a different location from <B|z|A> but they both have the same value.






Eigensystem[%%]
                   9     9   9  9
Out[8]= {{-9, 9, -(-), -(-), -, -, 0, 0, 0}, 
                   2     2   2  2
 
>    {{Sqrt[2], 0,  Sqrt[3], 0, 0, 0, 1, 0, 0}, 
 
>     {Sqrt[2], 0, -Sqrt[3], 0, 0, 0, 1, 0, 0}, 

>     {0, 0, 0, 1, 0, 0, 0, 1, 0}, 
 
>     {0, 1, 0, 0, 0, 1, 0, 0, 0}, 

>     {0, 0, 0, -1, 0, 0, 0, 1, 0}, 
 
>     {0, -1, 0, 0, 0, 1, 0, 0, 0}, 

>     {0, 0, 0, 0, 0, 0, 0, 0, 1}, 
            1
>     {-(-------), 0, 0, 0, 0, 0, 1, 0, 0},
         Sqrt[2]

>     {0, 0, 0, 0, 1, 0, 0, 0, 0}}}
