ContourPlot3DFor[i=0,i<11,i++,f[i,x_]=HermiteH[i,x] Exp[-x^2/2]/Sqrt[Sqrt[N[Pi]] 2^i i! ]]
f[16,x_]=HermiteH[16,x] Exp[-x^2/2]/Sqrt[Sqrt[N[Pi]] 2^16 16! ]
f[32,x_]=HermiteH[32,x] Exp[-x^2/2]/Sqrt[Sqrt[N[Pi]] 2^32 32! ]
Plot[{(f[0,x])^2+1,(f[1,x])^2+3,(f[2,x])^2+5,(f[3,x])^2+7,(f[4,x])^2+9,
x^2,1,3,5,7,9},{x,-4,4}, \
PlotRange->{0,11},PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
 {RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,1]},{RGBColor[1,0,0]},     
 {RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]

Plot[{(f[0,x])+1,(f[1,x])+3,(f[2,x])+5,(f[3,x])+7,(f[4,x])+9,
x^2,1,3,5,7,9},{x,-4,4}, \
PlotRange->{0,11},PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
 {RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,1]},{RGBColor[1,0,0]},     
 {RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]

Plot[{1,3,5,7,9,11,13,15,17,19},{x,0,1}, PlotRange->{-.05,20}]
Plot[f[32,x]^2,{x,-10.5,10.5},PlotPoints->150,AspectRatio->.67/GoldenRatio]
Plot[{f[32,x]^2,1/(Pi Sqrt[65-x^2]),{x,-10.5,10.5},PlotPoints->150,AspectRatio->.67/GoldenRatio]
Plot[1/(Pi Sqrt[65-x^2]),{x,-8,8},PlotStyle->{RGBColor[1,0,0]},AspectRatio->.67/GoldenRatio]

f[32,x_]=HermiteH[32,x] Exp[-x^2/2]/Sqrt[Sqrt[N[Pi]] 2^32 32! ]
-------
For[i=0,i<18,i++,
f[i,x_]=HermiteH[42+i,x] Exp[-x^2/2]/N[Sqrt[Sqrt[Pi] 2^(42+i) (42+i)! ]] ]

x0=x /. N[FindRoot[HermiteH[52,x],{x,9.4},AccuracyGoal->12],11]
{x -> 9.38527496322}
norm=1/Sqrt[NIntegrate[f[10,x]^2,{x,x0,13}]]
For[i=0,i<18,i++,b[i]=norm*NIntegrate[f[i,x] f[10,x],{x,x0,13}]]
bs=0
For[i=0,i<18,i++, Print[{i,b[i]}]; bs=bs+b[i]^2]
Print[bs]

p[t_,x_]=Sum[b[i]*f[i,x]*Exp[2 I(10-i)t],{i,0,17}]

Plot[Abs[p[0,x]]^2,{x,0,13},PlotRange->{0,1.5}]
Plot[Abs[p[.2,x]]^2,{x,0,13},PlotRange->{0,1.5}]
Plot[Abs[p[.4,x]]^2,{x,0,13},PlotRange->{0,.5}]
Plot[Abs[p[.6,x]]^2,{x,0,13},PlotRange->{0,.5}]
Plot[Abs[p[.8,x]]^2,{x,-6,6},PlotRange->{0,.5}]
Plot[Abs[p[1,x]]^2,{x,-13,0},PlotRange->{0,.5}]
Plot[Abs[p[1.2,x]]^2,{x,-13,0},PlotRange->{0,.5}]
Plot[Abs[p[1.4,x]]^2,{x,-13,0},PlotRange->{0,1.5}]
Plot[Abs[p[1.6,x]]^2,{x,-13,0},PlotRange->{0,1.5}]
Plot[Abs[p[1.8,x]]^2,{x,-13,0},PlotRange->{0,1.5}]
Plot[Abs[p[2,x]]^2,{x,-13,0},PlotRange->{0,.5}]
Plot[Abs[p[2.2,x]]^2,{x,-6,6},PlotRange->{0,.5}]
Plot[Abs[p[2.4,x]]^2,{x,-6,6},PlotRange->{0,.5}]
Plot[Abs[p[2.6,x]]^2,{x,0,13},PlotRange->{0,.5}]
Plot[Abs[p[2.8,x]]^2,{x,0,13},PlotRange->{0,.5}]
Plot[Abs[p[3,x]]^2,{x,0,13},PlotRange->{0,1.5}]
Plot[Abs[p[3.2,x]]^2,{x,0,13},PlotRange->{0,1.5}]

x0[0]=8.4
x1[0]=11.6
x0[1]=6
x1[1]=11.6
x0[2]=1
x1[2]=11.6
x0[3]=-2
x1[3]=10
x0[4]=-6
x1[4]=6
x0[5]=-10
x1[5]=1.5
x0[6]=-11.6
x1[6]=-2
x0[7]=-11.2
x1[7]=-6.8
x0[8]=-11.6
x1[8]=-8.4
x0[9]=-11.2
x1[9]=-6.2
x0[10]=-11.2
x1[10]=0
x0[11]=-10
x1[11]=4
x0[12]=-6
x1[12]=7
x0[13]=-1
x1[13]=10
x0[14]=3
x1[14]=11.2
x0[15]=7
x1[15]=11.2
x0[16]=8
x1[16]=11.6

For[i=0,i<17,i++,
pnorm[i]=NIntegrate[Abs[p[.2 i,x]]^2,{x,x0[i],x1[i]}]; Print[pnorm[i]] ]
For[i=0,i<17,i++,
xnorm[i]=NIntegrate[x Abs[p[.2 i,x]]^2,{x,x0[i],x1[i]}]; Print[xnorm[i]] ]
For[i=0,i<17,i++,
Print[{i,xnorm[i],pnorm[i],xnorm[i]/pnorm[i] }]]

res=Table[{.2 i,xnorm[i]/pnorm[i] },{i,0,16}]
Plot[10 Cos[2 x],{x,0,3.2},PlotStyle->{RGBColor[1,0,0]}]
ListPlot[res,PlotStyle->{AbsolutePointSize[5]}]
Show[%%,%]
------
dx0=NIntegrate[Abs[p[0,x]]^2 (x-xnorm[0]/pnorm[0])^2,{x,x0[0],x1[0]}]/pnorm[0]
0.0943743
p0[x_]=D[p[0,x],x]
dp0=NIntegrate[Abs[p0[x]]^2 ,{x,x0[0],x1[0]}]/pnorm[0]
3.66143

x00 =x /.{x -> 9.38527496322}
norm=1/Sqrt[NIntegrate[f[10,x]^2,{x,x00,13}]]
xavg=norm^2*NIntegrate[f[10,x] x f[10,x],{x,x00,13}]
xsqr=norm^2*NIntegrate[f[10,x] (x -xavg)^2 f[10,x],{x,x00,13}]
0.0650275
p00[x_]=D[f[10,x],x]
psqr=norm^2*NIntegrate[p00[x]^2,{x,x00,13}]
5.76995
Plot[norm f[10,x],{x,x00,13},PlotRange->{0,1.6}]
Plot[Abs[p[0,x]]^2,{x,8,13},PlotRange->{0,1.6},PlotStyle->{RGBColor[1,0,0]}]
Show[%%,%]
---
<<Statistics`ContinuousDistributions`
s[t_]=Sqrt[(.31 Cos[2 t])^2+(1.9 Sin[2 t])^2]
m[t_]=10 Cos[2 t]
pc[t_,x_]=PDF[NormalDistribution[m[t],s[t]],x]
Plot[{pc[0,x],pc[.2,x],pc[.4,x],pc[.6,x],pc[.8,x]},{x,-6,12},PlotRange->{0,1.5}]
---
p[n_,m_,r_]= r^(Abs[m]) LaguerreL[n,Abs[m],r^2] Exp[-r^2 /2] Sqrt[2/Pochhammer[n+1,Abs[m]]]
v[m_,r_]=m^2/r^2+r^2

Plot[{2p[0,3,x]+8,2p[1,3,x]+12,2p[2,3,x]+16,2p[3,3,x]+20,8,12,16,20,v[3,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2p[0,1,x]+4,2p[1,1,x]+8,2p[2,1,x]+12,2p[3,1,x]+16,2p[4,1,x]+20,4,8,12,16,20,v[1,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2p[0,0,x]+2,2p[1,0,x]+6,2p[2,0,x]+10,2p[3,0,x]+14,2p[4,0,x]+18,2p[5,0,x]+22,2,6,10,14,18,22,x^2},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

p1[i_,x_]=HermiteH[i,x] Exp[-x^2/2]/Sqrt[Sqrt[N[Pi]] 2^i i! ]

ContourPlot[p1[16,x]^2 p1[4,y]^2,{x,-8,8},{y,-4,4},PlotPoints->100,AspectRatio->.5]

ContourPlot[p[10,0,Sqrt[x^2+y^2]]^2,{x,-7,7},{y,-7,7},PlotPoints->100,ContourLines->False]
ContourPlot[p[5,10,Sqrt[x^2+y^2]]^2,{x,-7,7},{y,-7,7},PlotPoints->100,ContourLines->False]
---
Plot[{2,6,10,14},{x,-.15,.15}]
Plot[{6,10,14},{x,1.85,2.15}]
Plot[{6,10,14},{x,-2.15,-1.85}]
Plot[{10,14},{x,3.85,4.15}]
Plot[{10,14},{x,-4.15,-3.85}]
Plot[{14},{x,5.85,6.15}]
Plot[{14},{x,-6.15,-5.85}]
Plot[{4,8,12},{x,.85,1.15}]
Plot[{4,8,12},{x,-1.15,-.85}]
Plot[{8,12},{x,2.85,3.15}]
Plot[{8,12},{x,-3.15,-2.85}]
Plot[{12},{x,-5.15,-4.85}]
Plot[{12},{x,4.85,5.15}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,
%%%%%%%%%%,%%%%%%%%%%%,%%%%%%%%%%%%,%%%%%%%%%%%%%]
Show[%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%%,%%%%%%%%%%,
%%%%%%%%%%%%,%%%%%%%%%%%%,%%%%%%%%%%%%%,%%%%%%%%%%%%%%]
Plot[{2,4,6,8,10,12,14},{x,-.15,.15}]
Plot[{4,6,8,10,12,14},{x,.85,1.15}]
Plot[{6,8,10,12,14},{x,1.85,2.15}]
Plot[{8,10,12,14},{x,2.85,3.15}]
Plot[{10,12,14},{x,3.85,4.15}]
Plot[{12,14},{x,4.85,5.15}]
Plot[{14},{x,5.85,6.15}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%]
-----
LaguerreL[2,x^2+y^2]-(1/32) HermiteH[0,x] HermiteH[4,y]- 
(1/16) HermiteH[2,x] HermiteH[2,y] - (1/32) HermiteH[4,x] HermiteH[0,y] 

(x^2-y^2) LaguerreL[1,2,x^2+y^2]-(1/16) HermiteH[0,x] HermiteH[4,y]- 
(-1/16) HermiteH[4,x] HermiteH[0,y] 

(x^4-6x^2y^2+y^4)LaguerreL[0,4,x^2+y^2]-(1/16) HermiteH[0,x] HermiteH[4,y]- 
(-3/8) HermiteH[2,x] HermiteH[2,y] - (1/16) HermiteH[4,x] HermiteH[0,y] 

----
p[m_,r_,phi_]=r^m Exp[-r^2/2] Exp[I m phi] /Sqrt[Pi m!]
pm[r_,phi_]=N[Sqrt[1/(100 180!)]] r^180 Exp[-r^2/2]  \
Sum[ N[ Exp[-(m-20)^2/400]/Sqrt[Pochhammer[181,m]]] r^m Exp[ I m phi] ,{m,0,40}]
Plot[Abs[pm[14.1,phi]]^2,{phi,-Pi,Pi},PlotRange->{0,.21}]
ContourPlot[Abs[pm[Sqrt[x^2+y^2],ArcTan[x,y]]]^2,{x,12,16},{y,-2 ,2 },ContourLines->False]
ContourPlot[Abs[pm[Sqrt[x^2+y^2],ArcTan[x,y]]]^2,{x,-16,16},{y,-16,16},PlotPoints->100,ContourLines->False,
PlotRange->{0,.21}]
Plot[Abs[pm[x,0]]^2,{x,13,15.2},PlotRange->{0,.21}]
---
<<Graphics`ParametricPlot3D`
SphericalPlot3D[{1,GrayLevel[Abs[SphericalHarmonicY[5,0,theta,0]/SphericalHarmonicY[5,0,0,0]]^2]},
{theta, 0, Pi},{phi,0, 2 Pi},Lighting->False,PlotPoints->30]
Plot[Abs[SphericalHarmonicY[5,0,theta,0]]^2,{theta, 0, Pi},PlotPoints->50,PlotRange->{0,.9}]
SphericalPlot3D[{1,GrayLevel[Abs[SphericalHarmonicY[5,5,theta,0]/SphericalHarmonicY[5,5,Pi/2,0]]^2]},
{theta, 0, Pi},{phi,0, 2 Pi},Lighting->False,PlotPoints->30]
Plot[Abs[SphericalHarmonicY[5,5,theta,0]]^2,{theta, 0, Pi},PlotPoints->50,PlotRange->{0,.25}]
ParametricPlot3D[{Cos[p],Sin[p],0,RGBColor[1,0,0]},{p,0,2 Pi}]

Show[%,%%%,Lighting->False]
Show[%%,%%%%%%,Lighting->False]
----
p[n_,l_,m_,r_,theta_,phi_]=r^l Exp[-r^2/2] LaguerreL[n,l+1/2,r^2] SphericalHarmonicY[l,m,theta,phi] Sqrt[2 n!/(l+1/2+n)!]
pr[n_,l_,r_]=r^l Exp[-r^2/2] LaguerreL[n,l+1/2,r^2] Sqrt[2 n!/(l+1/2+n)!]
v[l_,x_]=l(l+1)/x^2+x^2

Plot[{2pr[0,3,x]+9,2pr[1,3,x]+13,2pr[2,3,x]+17,2pr[3,3,x]+21,9,13,17,21,v[3,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2pr[0,1,x]+5,2pr[1,1,x]+9,2pr[2,1,x]+13,2pr[3,1,x]+17,2pr[4,1,x]+21,5,9,13,17,21,v[1,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2pr[0,0,x]+3,2pr[1,0,x]+7,2pr[2,0,x]+11,2pr[3,0,x]+15,2pr[4,0,x]+19,3,7,11,15,19,x^2},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

ContourPlot[Abs[p[6,0,0,Sqrt[x^2+y^2],0,0]]^2,{x,-5,5},{y,-5,5},PlotPoints->100,ContourLines->False,
Contours->{.001,.002,.004,.008,.016,.032,.064,.128,.256,.512}]

Plot[Abs[p[6,0,0,x,0,0]]^2,{x,0,6},PlotPoints->50,PlotRange->{0,.55}]
Plot[50 Abs[p[6,0,0,x,0,0]]^2,{x,1.2,6},PlotPoints->50,PlotRange->{0,.55},PlotStyle->{RGBColor[0,0,1]}]
Show[%,%%]


Plot[Abs[p[2,2,2,r,Pi/2,0]]^2,{r,0,5},PlotPoints->50,PlotRange->{0,.075}]
ContourPlot[Abs[p[2,2,2,Sqrt[x^2+y^2],Pi/2,ArcTan[x,y]]]^2,{x,-5,5},{y,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.07}]
ContourPlot[Abs[p[2,2,2,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-5,5},{z,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.07}]

Plot[Abs[p[2,2,0,r,0,0]]^2,{r,0,5},PlotPoints->50,PlotRange->{0,.2}]
ContourPlot[Abs[p[2,2,0,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-5,5},{z,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.2}]

Plot[Abs[p[2,2,1,r,0,0]]^2,{r,0,5},PlotPoints->50,PlotRange->{0,.7}]
ContourPlot[Abs[p[2,2,1,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-5,5},{z,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.07}]

Plot[Abs[p[2,8,8,r,Pi/2,0]]^2,{r,0,6},PlotPoints->50,PlotRange->{0,.03}]
ContourPlot[Abs[p[2,8,8,Sqrt[x^2+y^2],Pi/2,ArcTan[x,y]]]^2,{x,-5,5},{y,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.03}]
ContourPlot[Abs[p[2,8,8,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-5,5},{z,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.03}]

Plot[Abs[p[2,8,0,r,0,0]]^2,{r,0,5},PlotPoints->50,PlotRange->{0,.15}]
ContourPlot[Abs[p[2,8,0,Sqrt[x^2+z^2],ArcTan[z,x],0]]^2,{x,-5,5},{z,-5,5},PlotPoints->100,ContourLines->False,
PlotRange->{0,.15}]

----
Plot[{3,7,11,15,19},{x,-.15,.15}]
Plot[{5,9,13,17},{x,.85,1.15}]
Plot[{7,11,15,19},{x,1.85,2.15}]
Plot[{9,13,17},{x,2.85,3.15}]
Plot[{11,15,19},{x,3.85,4.15}]
Plot[{13,17},{x,4.85,5.15}]
Plot[{15,19},{x,5.85,6.15}]
Plot[{17},{x,6.85,7.15}]
Plot[{19},{x,7.85,8.15}]
Show[%,%%,%%%,%%%%,%%%%%,%%%%%%,%%%%%%%,%%%%%%%%,%%%%%%%%,PlotRange->{0,20.5}]
-----
p[n_,l_,m_,r_,theta_,phi_]=r^l Exp[-r^2/2] LaguerreL[n,l+1/2,r^2] SphericalHarmonicY[l,m,theta,phi] Sqrt[2 n!/(l+1/2+n)!]
pr[n_,l_,r_]=r^(l+1) Exp[-r^2/2] LaguerreL[n,l+1/2,r^2] Sqrt[2 n!/(l+1/2+n)!]
v[l_,x_]=(l+.5)^2/x^2+x^2

Plot[{2pr[0,3,x]+9,2pr[1,3,x]+13,2pr[2,3,x]+17,2pr[3,3,x]+21,9,13,17,21,v[3,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2pr[0,1,x]+5,2pr[1,1,x]+9,2pr[2,1,x]+13,2pr[3,1,x]+17,2pr[4,1,x]+21,5,9,13,17,21,v[1,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[{2pr[0,0,x]+3,2pr[1,0,x]+7,2pr[2,0,x]+11,2pr[3,0,x]+15,2pr[4,0,x]+19,3,7,11,15,19,v[0,x]},{x,0,6},PlotRange->{0,25}, \
PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,0]},
{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,1]},{RGBColor[1,0,0]}}]

Plot[pr[32,0,x]^2,{x,0,13},PlotPoints->100]
Plot[(2/Pi)/Sqrt[ (131-.25/x^2-x^2)],{x,.05,11.4},PlotStyle->{RGBColor[1,0,0]}]
Show[%,%%,PlotRange->{0,.4}]
----
int1[l1_,m1_,l2_,m2_,l3_,m3_]:= (-1)^(m1)Integrate[Sin[theta]\
Integrate[ SphericalHarmonicY[l1,-m1,theta,phi]\
SphericalHarmonicY[l2,m2,theta,phi]\
SphericalHarmonicY[l3,m3,theta,phi],{phi,0,2 Pi}],{theta,0,Pi}]

int2[l1_,m1_,l2_,m2_,l3_,m3_]:=(-1)^(-m1)Sqrt[(2l1+1)(2l2+1)(2l3+1)/(4 Pi)]\
ThreeJSymbol[{l1,0},{l2,0},{l3,0}]ThreeJSymbol[{l1,-m1},{l2,m2},{l3,m3}]

int3[l1_,m1_,l2_,m2_,l3_,m3_]:=(-1)^(l_2-l_3)Sqrt[(2l_2+1)(2l_3+1)/(4 Pi)]\
ThreeJSymbol[{l1,0},{l2,0},{l3,0}]ClebschGordan[{l2,m2},{l3,m3},{l1,m1}]

int2[2,1,2,2,2,-1]
int2[2,2,2,2,2,0]
int2[2,2,2,2,0,0]

p[n_,l_,r_]=r^l Exp[-r^2/2] LaguerreL[n,l+1/2,r^2] Sqrt[2 n!/(l+1/2+n)!]

Integrate[p[1,0,r]^2 r^2 ,{r,0,Infinity}]

Integrate[p[1,0,r]^2 r^2 r^2 ,{r,0,Infinity}]
Integrate[p[1,0,r]p[0,2,r] r^2 r^2 ,{r,0,Infinity}]

veven={{0,-Sqrt[2/3], -Sqrt[4/3], 0},\
{-Sqrt[2/3],0,0,-Sqrt[2/3]},\
{-Sqrt[4/3],0,0,-Sqrt[4/3]},\
{0,-Sqrt[2/3], -Sqrt[4/3], 0}}
Eigensystem[veven]
% 
%{{-2, 0, 0, 2}, {{1, Sqrt[2/3], 2/Sqrt[3], 1}, {-1, 0, 0, 1},{0, -Sqrt[2], 1, 0}, {1, -Sqrt[2/3], -2/Sqrt[3], 1}}}
%
vodd={{0,-1},{-1,0}}
Eigensystem[vodd]
(SphericalHarmonicY[2,1,theta,phi]- SphericalHarmonicY[2,-1,theta,phi])/Sqrt[2]
ExpToTrig[%] (* Math 3 untested *)

%result is -Sqrt[15/Pi]Sin[2 theta]Cos[phi]/4

yB[theta_,phi_]:=-Sqrt[15/Pi]Sin[2 theta]Cos[phi]/4

Integrate[Sin[theta]\
Integrate[ yB[theta,phi]^2 (Sin[theta]^2 Cos[2 phi])^2,{phi,0,2 Pi}],{theta,0,Pi}]
Integrate[Sin[theta]\
Integrate[ yB[theta,phi]^2 (Sin[theta]^2 Cos[2 phi]),{phi,0,2 Pi}],{theta,0,Pi}]
Integrate[p[0,2,r]^2 r^4 r^2 ,{r,0,Infinity}]