Coordenadas rectangulares

La ecuación de Schrödinger independiente del tiempo para una partícula de masa m_e en el potencial armónico V(x,y)

${\displaystyle -\frac{{\hslash }^2}{2m_e}\left(\frac{{\partial }^2\psi (x,y)}{\partial x^2}+\frac{{\partial }^2\psi (x,y)}{\partial y^2}\right)+\frac{1}{2}{m}_e{\omega }^2\left(x^2+y^2\right)\psi (x,y)=E\psi (x,y)}$

Tomamos una escala de energías y distancias

${\displaystyle x=\xi \sqrt{\frac{\hslash }{m\omega }},y=\eta \sqrt{\frac{\hslash }{m\omega }},E=\hslash \omega \epsilon }$

La ecuación de Schrödinger se transforma en otra más simple

${\displaystyle \begin{array}{l}-\frac{{\partial }^2\psi (\xi ,\eta )}{\partial {\xi }^2}-\frac{{\partial }^2\psi (\xi ,\eta )}{\partial {\eta }^2}+\left({\xi }^2+{\eta }^2\right)\psi (\xi ,\eta )=2\epsilon \psi (\xi ,\eta )\\ \frac{{\partial }^2\psi (\xi ,\eta )}{\partial {\xi }^2}+\frac{{\partial }^2\psi (\xi ,\eta )}{\partial {\eta }^2}+\left\{2\epsilon -\left({\xi }^2+{\eta }^2\right)\right\}\psi (\xi ,\eta )=0\end{array}}$

Probamos la solución Ψ(ξ,η)=X(ξ)Y(η), separación de variables

${\displaystyle \begin{array}{l}Y(\eta )\frac{d^{2}X(\xi )}{d{\xi }^2}+X(\xi )\frac{d^{2}Y(\eta )}{d{\eta }^2}+\left\{2\epsilon -\left({\xi }^2+{\eta }^2\right)\right\}X(\xi )Y(\eta )=0\\ \frac{1}{X(\xi )}\frac{d^{2}X(\xi )}{d{\xi }^2}+\left\{2{\epsilon }_x-{\xi }^2\right\}+\frac{1}{Y(\eta )}\frac{d^{2}Y(\eta )}{d{\eta }^2}+\left\{2{\epsilon }_y-{\eta }^2\right\}=\mathrm{0,}{\epsilon }_x+{\epsilon }_y=\epsilon \end{array}}$

El primer término, es una función solamente de ξ, el segundo, de η. Esta ecuación diferencial se convierte en el sistema de dos ecuaciones diferenciales

${\displaystyle \{\begin{array}{l}\frac{d^{2}X(\xi )}{d{\xi }^2}+\left\{2{\epsilon }_x-{\xi }^2\right\}X(\xi )=0\\ \frac{d^{2}Y(\eta )}{d{\eta }^2}+\left\{2{\epsilon }_y-{\eta }^2\right\}Y(\eta )=0\end{array}}$

Haciendo un cambio de variable, la primera ecuación diferencial se transforma en la de Hermite

${\displaystyle \begin{array}{l}X(\xi )=x(\xi )\exp \left(-\frac{{\xi }^2}{2}\right)\hfill \\ \frac{dX}{d\xi }=\exp \left(-\frac{{\xi }^2}{2}\right)\frac{dx}{d\xi }-\xi \exp \exp \left(-\frac{{\xi }^2}{2}\right)x\hfill \\ \frac{d^{2}X}{d{\xi }^2}=\exp \left(-\frac{{\xi }^2}{2}\right)\frac{d^{2}x}{d{\xi }^2}-\xi \exp \left(-\frac{{\xi }^2}{2}\right)\frac{dx}{d\xi }-\exp \left(-\frac{{\xi }^2}{2}\right)x+{\xi }^2\exp \left(-\frac{{\xi }^2}{2}\right)x-\xi \exp \left(-\frac{{\xi }^2}{2}\right)\frac{dx}{d\xi }\hfill \\ \frac{d^{2}x}{d{\xi }^2}-2\xi \frac{dx}{d\xi }+(2{\epsilon }_x-1)x=0\hfill \end{array}}$

Del mismo modo, la segunda

${\displaystyle \frac{d^{2}y}{d{\eta }^2}-2\eta \frac{dy}{d\eta }+(2{\epsilon }_y-1)y=0}$

Cuya solución son los polinomios de Hermite

Los niveles de energía ε son

${\displaystyle \begin{array}{l}\{\begin{array}{l}2{\epsilon }_x-1=2n_x,{\epsilon }_x=n_x+\frac{1}{2},n_x=\mathrm{0,1,2,3...}\\ 2{\epsilon }_y-1=2n_y,{\epsilon }_y=n_y+\frac{1}{2},n_y=\mathrm{0,1,2,3...}\end{array}\\ {\epsilon }_{n_x,n_y}={\epsilon }_x+{\epsilon }_x=n_x+n_y+1\end{array}}$

Cuando n_x=1 y n_y=0, la energía es ε_1,0=2. Cuando n_x=0 y n_y=1, la energía es la misma. Dos estados con la misma energía

Las funciones de onda

${\displaystyle \begin{array}{l}x(\xi )=C·H_{n_x}(\xi ),y(\eta )=C·H_{n_y}(\eta )\\ {\psi }_{n_x,n_y}\left(\xi ,\eta \right)=C·H_{n_x}(\xi )H_{n_y}(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)\end{array}}$

Donde C es una constante a determinar

${\displaystyle \begin{array}{l}{C}^2{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{H}_{n_x}^2(\xi )\exp \left(-{\xi }^2\right)d\xi }·{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{H}_{n_y}^2(\eta )\exp \left(-{\eta }^2\right)d\eta =1}\hfill \\ C^2\left(\sqrt{\pi }{2}^{n_x}{n}_x!\right)\left(\sqrt{\pi }{2}^{n_y}{n}_y!\right)=1\hfill \end{array}}$

Se ha tenido en cuenta las relaciones de ortogonalidad de los polinomios de Hermite. El resultado final es

${\displaystyle {\psi }_{n_x,n_y}\left(\xi ,\eta \right)=\frac{1}{\sqrt{\pi ·2^{n_x+n_y}{n}_x!n_y!}}{H}_{n_x}(\xi )H_{n_y}(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)}$

Comprobaremos en los ejemplos que

${\displaystyle {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\left|{\psi }_{n_x,n_y}\left(\xi ,\eta \right)\right|}^{2}d\xi ·d\eta =\frac{1}{\pi ·2^{n_x+n_y}{n}_x!n_y!}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{H}_{n_x}^2(\xi )}}\exp \left(-{\xi }^2\right)d\xi ·{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{H}_{n_x}^2(\xi )\exp \left(-{\eta }^2\right)d\eta }}=1}$

Ejemplos

• Ψ_0,0(ξ,η)

${\displaystyle {\psi }_{\mathrm{0,0}}\left(\xi ,\eta \right)=\frac{1}{\sqrt{\pi ·2^0·0!·0!}}{H}_0(\xi )H_0(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)=\frac{1}{\sqrt{\pi }}\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)}$

[x,y] = meshgrid(-5:0.2:5);
z = exp(-(x.^2+y.^2)/2)/sqrt(pi);
surf(x,y,z)
xlabel('x'); 
ylabel('y'); 
zlabel('z')
title('\Psi_{0,0}')



Comprobación.

>> syms x y;
>> int(exp(-x^2),x,-inf,inf)*int(exp(-y^2),y,-inf,inf)/pi
 ans =1

• Ψ_1,0(ξ,η)

${\displaystyle {\psi }_{\mathrm{1,0}}\left(\xi ,\eta \right)=\frac{1}{\sqrt{\pi ·2^1·1!·0!}}{H}_1(\xi )H_0(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)=\frac{1}{\sqrt{2\pi }}2\xi \exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)}$

[x,y] = meshgrid(-5:0.2:5);
z = x.*exp(-(x.^2+y.^2)/2)*sqrt(2/pi);
surf(x,y,z)
xlabel('x'); 
ylabel('y'); 
zlabel('z')
title('\Psi_{1,0}')



Comprobación.

>> syms x y;
>> int(x^2*exp(-x^2),x,-inf,inf)*int(exp(-y^2),y,-inf,inf)*2/pi
 ans =1

• Ψ_2,0(ξ,η)

${\displaystyle {\psi }_{\mathrm{2,0}}\left(\xi ,\eta \right)=\frac{1}{\sqrt{\pi ·2^2·2!·0!}}{H}_2(\xi )H_0(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)=\frac{1}{2\sqrt{2\pi }}\left(4{\xi }^2-2\right)\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)}$

[x,y] = meshgrid(-5:0.2:5);
z = (4*x.^2-2).*exp(-(x.^2+y.^2)/2)/sqrt(8*pi);
surf(x,y,z)
xlabel('x'); 
ylabel('y'); 
zlabel('z')
title('\Psi_{2,0}')



Comprobación.

>> syms x y;
>> int((4*x^2-2)^2*exp(-x^2),x,-inf,inf)*int(exp(-y^2),y,-inf,inf)/(8*pi)
 ans =1

• Ψ_2,1(ξ,η)

${\displaystyle {\psi }_{\mathrm{2,1}}\left(\xi ,\eta \right)=\frac{1}{\sqrt{\pi ·2^3·2!·1!}}{H}_2(\xi )H_1(\eta )\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)=\frac{1}{2\sqrt{\pi }}\eta \left(4{\xi }^2-2\right)\exp \left(-\frac{{\xi }^2+{\eta }^2}{2}\right)}$

[x,y] = meshgrid(-5:0.2:5);
z = y.*(4*x.^2-2).*exp(-(x.^2+y.^2)/2)/sqrt(4*pi);
surf(x,y,z)
xlabel('x'); 
ylabel('y'); 
zlabel('z')
title('\Psi_{2,1}')



Comprobación.

>> syms x y;
>> int((4*x^2-2)^2*exp(-x^2),x,-inf,inf)*int(y^2*exp(-y^2),y,-inf,inf)/(4*pi)
 ans =1

Coordenadas polares

Resolveremos la ecuación de Schrödinger independiente del tiempo en coordenadas polares, para el potencial V(r). Donde m_e es la masa de la partícula

${\displaystyle -\frac{{\hslash }^2}{2m_e}\left(\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \psi (r,\phi )}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\psi (r,\phi )}{\partial {\phi }^2}\right)+\left(\frac{1}{2}{m}_e{\omega }^2{r}^2\right)\psi (r,\phi )=E\psi (r,\phi )}$

Intentamos la separación de variables

${\displaystyle \begin{array}{l}\psi (r,\phi )=R(r)·\Phi (\phi )\\ \frac{\Phi }{r}\frac{d}{dr}\left(r\frac{dR}{dr}\right)+\frac{R}{r^2}\frac{d^2\Phi }{d{\phi }^2}+\frac{2m_e}{{\hslash }^2}\left(E-\frac{1}{2}m{\omega }^2{r}^2\right)R\Phi =0\\ \frac{r}{R}\frac{d}{\partial r}\left(r\frac{dR}{dr}\right)+\frac{2m_e}{{\hslash }^2}\left(E-\frac{1}{2}{m}_e{\omega }^2{r}^2\right)r^2=-\frac{1}{\Phi }\frac{d^2\Phi }{d{\phi }^2}\end{array}}$

El primer término, solamente depende de r, el segundo de φ, la ecuación diferencial se transforma en el sistema de dos ecuaciones diferenciales

${\displaystyle \begin{array}{l}\{\begin{array}{l}\frac{r}{R}\frac{d}{dr}\left(r\frac{dR}{dr}\right)+\frac{2m_e}{{\hslash }^2}\left(E-\frac{1}{2}{m}_e{\omega }^2{r}^2\right)r^2=m^2\\ -\frac{1}{\Phi }\frac{d^2\Phi }{d{\phi }^2}=m^2\end{array}\\ \{\begin{array}{l}\frac{d^{2}R}{d{r}^2}+\frac{1}{r}\frac{dR}{dr}+\left\{\frac{2m_e}{{\hslash }^2}\left(E-\frac{1}{2}{m}_e{\omega }^2{r}^2\right)-\frac{m^2}{r^2}\right\}R=0\\ \frac{d^2\Phi }{d{\phi }^2}+m^2\Phi =0\end{array}\end{array}}$

La ecuación angular

Tiene una solución conocida, que expresamos de forma equivalente

${\displaystyle \Phi (\phi )=C_1\exp (im\phi )+C_2\exp (-im\phi )}$

La solución es periódica, Φ(φ+2π)=Φ(φ). m tiene que ser un entero.

${\displaystyle \Phi (\phi )=C\exp (im\phi ),m=\mathrm{0,}\pm \mathrm{1,}\pm \mathrm{2....}}$

La ecuación radial

Haciendo el cambio de variable, la ecuación radial se transforma en

${\displaystyle \begin{array}{l}E=\hslash \omega \epsilon ,r=\sqrt{\frac{\hslash }{m_e\omega }}x\\ \frac{d^{2}R}{d{x}^2}+\frac{1}{x}\frac{dR}{dx}+\left\{2\epsilon -x^2-\frac{m^2}{x^2}\right\}R=0\end{array}}$

Definimos

${\displaystyle R(x)=x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right)G(x)}$

Calculamos la derivada primera y segunda de R respecto de x

${\displaystyle \begin{array}{l}\frac{dR}{dx}=\left|m\right|x^{\left|m\right|-1}\exp \left(-\frac{x^2}{2}\right)G(x)-x^{\left|m\right|+1}\exp \left(-\frac{x^2}{2}\right)G(x)+x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right)\frac{dG(x)}{dx}\\ \frac{d^{2}R}{d{x}^2}=\left|m\right|\left(\left|m\right|-1\right)x^{\left|m\right|-2}\exp \left(-\frac{x^2}{2}\right)G(x)-\left|m\right|x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right)G(x)+\left|m\right|x^{\left|m\right|-1}\exp \left(-\frac{x^2}{2}\right)\frac{dG(x)}{dx}\\ -\left(\left|m\right|+1\right)x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right)G(x)+x^{\left|m\right|+2}\exp \left(-\frac{x^2}{2}\right)G(x)-x^{\left|m\right|+1}\exp \left(-\frac{x^2}{2}\right)\frac{dG(x)}{dx}\\ +\left|m\right|x^{\left|m\right|-1}\exp \left(-\frac{x^2}{2}\right)\frac{dG(x)}{dx}-x^{\left|m\right|+1}\exp \left(-\frac{x^2}{2}\right)\frac{dG(x)}{dx}+x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right)\frac{d^{2}G(x)}{d{x}^2}\end{array}}$

Obtenemos una ecuación diferencial en G

${\displaystyle \begin{array}{l}{x}^{\left|m\right|}\frac{d^{2}G(x)}{d{x}^2}+\left\{\left(2\left|m\right|+1\right)x^{\left|m\right|-1}-2x^{\left|m\right|+1}\right\}\frac{dG(x)}{dx}+\left\{2\epsilon -2\left(\left|m\right|+1\right)\right\}{x}^{\left|m\right|}G(x)=0\\ \frac{d^{2}G(x)}{d{x}^2}+\left(\frac{\left(2\left|m\right|+1\right)}{x}-2x\right)\frac{dG(x)}{dx}+\left(2\epsilon -2\left(\left|m\right|+1\right)\right)G(x)=0\end{array}}$

Hacemos el cambio de variable ξ=x²

${\displaystyle \begin{array}{l}\frac{dG}{dx}=\frac{dG}{d\xi }\frac{d\xi }{dx}=2x\frac{dG}{d\xi }\\ \frac{d^{2}G}{d{x}^2}=\frac{d}{dx}\left(2x\frac{dG}{d\xi }\right)=2\frac{dG}{d\xi }+2x\frac{d}{dx}\left(\frac{dG}{d\xi }\right)=2\frac{dG}{d\xi }+2x\frac{d}{d\xi }\left(\frac{dG}{d\xi }\right)\frac{d\xi }{dx}=2\frac{dG}{d\xi }+4x^2\frac{d^{2}G}{d{\xi }^2}\\ 2\frac{dG}{d\xi }+4x^2\frac{d^{2}G}{d{\xi }^2}+\left(\frac{2\left|m\right|+1}{x}-2x\right)2x\frac{dG}{d\xi }+\left(2\epsilon -2\left(\left|m\right|+1\right)\right)G=0\\ 2\frac{dG}{d\xi }+4\xi \frac{d^{2}G}{d{\xi }^2}+\left(4\left|m\right|+2-4\xi \right)\frac{dG}{d\xi }+\left(2\epsilon -2\left(\left|m\right|+1\right)\right)G=0\\ \xi \frac{d^{2}G(\xi )}{d{\xi }^2}+\left(\left|m\right|+1-\xi \right)\frac{dG(\xi )}{d\xi }-\frac{\left|m\right|+1-\epsilon }{2}G(\xi )=0\end{array}}$

Es una ecuación del tipo

${\displaystyle z\frac{d^{2}w}{d{z}^2}+\left(b-z\right)\frac{dw}{dz}-aw=0}$

que se denomina hipergeométrica confluente, cuya solución es

${\displaystyle \begin{array}{l}G(\xi )={}_{1}F{}_1\left(a;b,\xi \right)={}_{1}F{}_1\left(-n;\left|m\right|+\mathrm{1,}\xi \right),-n=\frac{\left|m\right|+1-\epsilon }{2}\\ G(x^2)={}_{1}F{}_1\left(-n;\left|m\right|+\mathrm{1,}{x}^2\right)\\ R_{n,m}(x)=C·x^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-n;\left|m\right|+\mathrm{1,}{x}^2\right)\end{array}}$

n es un número entero positivo (incluido cero) que se denomina número cuántico principal. Los niveles de energía se han obtenido en el primer apartado, ε=n_x+n_y+1=k+1, k=0,1,2,3...

${\displaystyle \begin{array}{l}-n=\frac{\left|m\right|+1-\epsilon }{2},k=\mathrm{0,1,2,3...}\\ -n=\frac{\left|m\right|-k}{2},n=\mathrm{0,1,2,3...,}k=\mathrm{0,1,2,3...}\end{array}}$

Dado el entero k, para que n sea entero, los posibles valores de |m| se recogen en la tabla

Las funciones de onda radiales R_nm(r) son

${\displaystyle R_{nm}(x)=\frac{1}{\left|m\right|!}\sqrt{\frac{2\left(n+\left|m\right|\right)!}{n!}}{x}^{\left|m\right|}\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-n;\left|m\right|+\mathrm{1,}{x}^2\right)}$

Esta función hipergeométrica está relacionada con los polinomios asociados de Laguerre del siguiente modo

${\displaystyle L_n^m(z)=\frac{\left(n+m\right)!}{n!m!}{}_{1}F{}_1\left(-n;m+\mathrm{1,}z\right)}$

Funciones radiales R_nm(r)

Comprobaremos utilizando la función int de MATLAB que las funciones radiales están normalizadas

La integral

${\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}{R}_{nm}^2(r)\left(2\pi r·dr\right)}=1}$

también

${\displaystyle {\displaystyle \underset{0}{\overset{\infty }{\int }}r{R}_{nm}^2(r)·dr}=1}$

• R₀₀(r)

${\displaystyle \begin{array}{l}{R}_{00}(x)=\frac{1}{0!}\sqrt{\frac{2\left(0+0\right)!}{0!}}{x}^0\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-1;0+\mathrm{1,}{x}^2\right)\\ R_{00}(r)=\sqrt{2}\exp \left(-\frac{x^2}{2}\right)\frac{0!0!}{\left(0+0\right)!}{L}_0^0(x^2)=\sqrt{2}\exp \left(-\frac{x^2}{2}\right)\end{array}}$

Representamos R₀₀(x) en función de x

f=@(x) sqrt(2)*exp(-x.^2/2);
fplot(f,[0,5])
grid on
xlabel('x')
ylabel('R_{00}(x)')
title('Funciones radiales, R_{0,0}(x)')



>> syms x;
>> int(x*exp(-x^2),x,0,inf)*2
 ans =1

• R₁₁(x)

${\displaystyle \begin{array}{l}{R}_{11}(x)=\frac{1}{1!}\sqrt{\frac{2\left(1+1\right)!}{1!}}{x}^1\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-1;1+\mathrm{1,}{x}^2\right)\\ R_{11}(r)=2x\exp \left(-\frac{x^2}{2}\right)\frac{1!1!}{\left(1+1\right)!}{L}_1^1(x^2)=x\exp \left(-\frac{x^2}{2}\right)\left(2-x^2\right)\end{array}}$

f=@(x) (x.*exp(-x.^2/2)).*(2-x.^2);
fplot(f,[0,5])
grid on
xlabel('x')
ylabel('R_{11}(x)')
title('Funciones radiales, R_{1,1}(x)')



>> syms x;
>> int(x^3*(2-x^2)^2*exp(-x^2),x,0,inf)
 ans =1

• R₂₀(x)

${\displaystyle \begin{array}{l}{R}_{20}(x)=\frac{1}{0!}\sqrt{\frac{2\left(2+0\right)!}{2!}}{x}^0\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-2;0+\mathrm{1,}{x}^2\right)\\ R_{20}(r)=\sqrt{2}\exp \left(-\frac{x^2}{2}\right)\frac{2!0!}{\left(2+0\right)!}{L}_2^0(x^2)=\sqrt{2}\exp \left(-\frac{x^2}{2}\right)\left(\frac{1}{2}{x}^4-2x^2+1\right)\end{array}}$

• R₂₂(x)

${\displaystyle \begin{array}{l}{R}_{22}(x)=\frac{1}{2!}\sqrt{\frac{2\left(2+2\right)!}{2!}}{x}^2\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-2;2+\mathrm{1,}{x}^2\right)\\ R_{22}(x)=\sqrt{6}\exp \left(-\frac{x^2}{2}\right)x^2\frac{2!·2!}{\left(2+2\right)!}{L}_2^2(x^2)=\frac{1}{\sqrt{6}}\exp \left(-\frac{x^2}{2}\right)x^2\left(\frac{1}{2}{x}^4-4x^2+6\right)\end{array}}$

hold on
f=@(x) sqrt(2)*exp(-x.^2/2).*(x.^4/2-2*x.^2+1);
fplot(f,[0,5])
g=@(x) (x.^2.*exp(-x.^2/2)).*(x.^4/2-4*x.^2+6)/sqrt(6);
fplot(g,[0,5])
hold off
grid on
xlabel('x')
ylabel('R_{20}(x), R_{2,2}(x)')
title('Funciones radiales, R_{2,0}(x), R_{2,2}(x)')



>> syms x;
>> 2*int(x*(1-x^2)^2*exp(-x^2),x,0,inf)
 ans =1
>> int(x^5*(x^4/2-4*x^2+6)^2*exp(-x^2),x,0,inf)/6
 ans =1

• R₃₁(x)

${\displaystyle \begin{array}{l}{R}_{31}(x)=\frac{1}{1!}\sqrt{\frac{2\left(3+1\right)!}{3!}}{x}^1\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-3;1+\mathrm{1,}{x}^2\right)\\ R_{31}(x)=2\sqrt{2}\exp \left(-\frac{x^2}{2}\right)x\frac{3!·1!}{\left(3+1\right)!}{L}_3^1(x^2)=\frac{\sqrt{2}}{12}\exp \left(-\frac{x^2}{2}\right)x\left(-x^6+12x^4-36x^2+24\right)\end{array}}$

• R₃₃(x)

${\displaystyle \begin{array}{l}{R}_{33}(x)=\frac{1}{3!}\sqrt{\frac{2\left(3+3\right)!}{3!}}{x}^3\exp \left(-\frac{x^2}{2}\right){}_{1}F{}_1\left(-3;3+\mathrm{1,}{x}^2\right)\\ R_{33}(x)=\sqrt{\frac{20}{3}}\exp \left(-\frac{x^2}{2}\right)x^3\frac{3!·3!}{\left(3+3\right)!}{L}_3^3(x^2)=\frac{1}{6\sqrt{60}}\exp \left(-\frac{x^2}{2}\right)x^3\left(-x^6+18x^4-90x^2+120\right)\end{array}}$

hold on
f1=@(x) sqrt(2)*(exp(-x.^2/2).*x).*(-x.^6+12*x.^4-36*x.^2+24)/12;
fplot(f1,[0,5])
f2=@(x) (exp(-x.^2/2).*x.^3).*(-x.^6+18*x.^4-90*x.^2+120)/(6*sqrt(60));
fplot(f2,[0,5])
hold off
grid on
xlabel('x')
ylabel('R_{3,1}(x), R_{3,3}(x)')
title('Funciones radiales, R_{3,1}(x), R_{3,3}(x)')

>> syms x;
>> int(x^3*(-x^6+12*x^4-36*x^2+24)^2*exp(-x^2),x,0,inf)/72
 ans =1
>> int(x^7*(-x^6+18*x^4-90*x^2+120)^2*exp(-x^2),x,0,inf)/(36*60)
 ans =1

Alternativamente, utilizamos la función laguerreL(n,m,z) de MATLAB en vez de los polinomios asociados de Laguerre $L_n^m(z)$

hold on
f1=@(x) sqrt(2)*(exp(-x.^2/2).*x).*laguerreL(3,1,x.^2)/2;
fplot(f1,[0,5])
f2=@(x) (exp(-x.^2/2).*x.^3).*laguerreL(3,3,x.^2)/sqrt(60);
fplot(f2,[0,5])
hold off
grid on
xlabel('x')
ylabel('R_{3,1}(x), R_{3,3}(x)')
title('Funciones radiales, R_{3,1}(x), R_{3,3}(x)')

>> int((exp(-x.^2).*x.^3).*laguerreL(3,1,x.^2).^2/2,x,0,inf)
 ans =1
 >> int((exp(-x.^2).*x.^7).*laguerreL(3,3,x.^2).^2/60,x,0,inf)
 ans =1

Referencias

Xijin Fu. Study on two - dimensional linear harmonic oscillator characteristics based on MATLAB software. IOP Conf. Series: Earth and Environmental Science 295(2019) 032042. IOP Publishing