Атом водорода.
Собственные состояния атома водорода задаются числами l,m,n
n-главное квантовое число, n - натуральное
l-квантовое число углового момента может принимать значения 0,1,...,n-1  и определяет величину углового момента
m-магнитное  квантовое число может примать значения-l,-l+1,...,0,...l и определяет проекцию момента импульса на ось z
Y[m,l,θ, ϕ] -сферические функции
Xi[r,n, l] - радиальные волновые функции
Собственные состояния имеют вид
|m,l,n>= 1/r Xi[r,n, l]Y[m,l,θ, ϕ] c энергией E_n=-1/(2(n^2))m e^4
Уровни с данной энергией E_n      (n^2)-кратно вырождены.
Квадрат модуля волновых функций определяет плотность вероятности обнаружить частицу в элементе объема r^2 Sin[θ]drdθdϕ.
Далее Probden1[θ, ϕ] - плотность вероятности обнаружить частицу в элементе телесного угла Sin[θ}dθdϕ,
Probden2[u] - плотность вероятности обнаружить частицу в интервале du.

Clear[l,m,θ, ϕ]
Z[0,l_,θ_, ϕ_]:=Sin[θ]^l Exp[I l ϕ];
Z[a_,l_,θ_, ϕ_]:= -Exp[-I ϕ]D[Z[a-1,l,θ, ϕ],θ]+I Exp[-I ϕ]Cot[θ]D[Z[a-1,l,θ, ϕ],ϕ];
Con[m_,l_]:=Sqrt[((2l+1)(l+m)!)/(π (l-m)!)]1/(l!2^(l+1));
Y[m_,l_,θ_, ϕ_]:=Simplify[Con[m,l] Z[l-m,l,θ, ϕ]];
Ym2[m_,l_,θ_, ϕ_]:=(Exp[-m I ϕ]Y[m,l,θ, ϕ])^2;
Probden1[θ_, ϕ_]=Ym2[0,1,θ, ϕ]
SphericalPlot3D[Probden1[θ, ϕ],{θ,0,π},{ ϕ,0,2π}]

SphericalPlot3D[((5 (((1+3 Cos[2 θ]))^2))/(64 π)),{θ,0,π},{ ϕ,0,2π}]
SphericalPlot3D[((15 ((Cos[θ])^2) ((Sin[θ])^2))/(8 π)),{θ,0,π},{ ϕ,0,2π}]
SphericalPlot3D[((15 ((Sin[θ])^4))/(32 π)),{θ,0,π},{ ϕ,0,2π}]
(*Это плотность вероятности обнаружить частицу в элементе телесного угла Sin[θ}dθdϕ для d-орбиталей*)

Clear[n,l];
cc[0,n_,l_]:=1;
cc[k_,n_,l_]:=2*(k+l-n)*cc[k-1,n,l]/(k-1+2l+2)/k/n;
Lag[u_,n_,l_]:=Sum[cc[i,n,l] u^i,{i,0,n-l-1}];
Nn[n_,l_]:=Integrate[u^(2l+2)Exp[-2u/n]Lag[u,n,l]^2,{u,0,Infinity}];
Xi[u_,n_, l_]:=1/Simplify[Sqrt[Nn[n,l]]]u^(l+1)Exp[-u/n]Lag[u,n,l];
Probden2[u_]=Xi[u,4,1]^2
Plot[Probden2[u],{u,0,50}];
(*Например, плотность вероятности для 4p-орбитали*)

Plot[{(4/27) (E^(-2 u/3)) (u^2) (((1-((2 u)/3)+((2 (u^2))/27)))^2),((32 (E^(-2 u/3)) (((1-(u/6)))^2) (u^4))/2187),((8 (E^(-2 u/3)) (u^6))/98415)},{u,0,50}]
(*На рисунке изображена плотность вероятности обнаружить частицу в интервале du для 3s,3p,3d орбиталей *)

Clear[m,l,n]
EigFun[m_,l_,n_,θ_, ϕ_,r_]:=1/r Xi[r,n, l]Y[m,l,θ, ϕ];
EigFun[1,2,3,θ, ϕ,r];
m=0;
l=0;
n=4;
Probden[θ_, ϕ_,r_]=(Exp[-m I  ϕ]EigFun[m,l,n,θ,  ϕ,r])^2
Plot3D[Probden[ArcTan[y/x],π,Sqrt[x^2+y^2]],{x,0.1,50},{y,-80,80}]
(*Это плотность вероятности на срезе, то есть при фиксированном угле ϕ*)