I have a question regarding simulating the hydrogen atom and the atomic orbitals.

So I wanted to mess around a little bit with sliders in Python and I figured, that an animation of the atomic orbitals for different n, l and m values would be ideal. The part with the sliders was easy, and I could just stop there, but I am too determined in figuring out why I get wrong results.
In the following, you can find my code:

import numpy as np
from scipy.special import sph_harm, eval_genlaguerre
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider, Button
import math

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection='3d')

# Initial values for the quantum numbers
n_init = 1
l_init = 0
m_init = 0

# Arrays for the spherical coordinates
phi = np.linspace(0, 2 * np.pi, 150)
theta = np.linspace(0, np.pi, 150)
R = np.linspace(0, 20, 150)

PHI, THETA = np.meshgrid(phi, theta)

surf = ax.plot_surface(np.array([[]]), np.array([[]]), np.array([[]]), color='0.75')

ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)

# Sliders for the quantum numbers
n_freq = fig.add_axes([0.25, 0.15, 0.65, 0.03])
l_freq = fig.add_axes([0.25, 0.1, 0.65, 0.03])
m_freq = fig.add_axes([0.25, 0.05, 0.65, 0.03])
n_slider = Slider(ax=n_freq, label='n', valmin=1, valmax=10, valinit=n_init, valstep=1)
l_slider = Slider(ax=l_freq, label='l', valmin=0, valmax=n_slider.valmax - 1, valinit=l_init, valstep=1)
m_slider = Slider(ax=m_freq, label='m', valmin=-l_slider.valmax, valmax=l_slider.valmax, valinit=m_init, valstep=1)

# Update the plot when the sliders are changed
def update(val):
    n = int(n_slider.val)
    l = int(l_slider.val)
    
    l_slider.valmax = n - 1
    if l > l_slider.valmax:
        l_slider.set_val(l_slider.valmax)
        l = l_slider.valmax
        
    m_slider.valmax = l
    m_slider.valmin = -l
    if int(m_slider.val) > l:
        m_slider.set_val(l)
    if int(m_slider.val) < -l:
        m_slider.set_val(-l)

    fig.canvas.draw_idle()

n_slider.on_changed(update)
l_slider.on_changed(update)
m_slider.on_changed(update)

# Introducing a reset button for the sliders
resetax = fig.add_axes([0.08, 0.065, 0.1, 0.04])
button = Button(resetax, 'Reset', hovercolor='0.975')

def reset(event):
    l_slider.reset()
    m_slider.reset()
    n_slider.reset()
    
button.on_clicked(reset)

# Function to calculate the radial component of the wave function
def R_nl(r, n, l):
    a_0 = 1
    rho = 2 * r / (n * a_0)
    normalization = np.sqrt((2 / (n * a_0))**3 * math.factorial(n - l - 1) / (2 * n * math.factorial(n + l)))
    radial_component = normalization * np.exp(-rho / 2) * rho**l * eval_genlaguerre(n - l - 1, 2 * l + 1, rho)
    return radial_component

def animate(i):
    global surf
    t = 2 * np.pi * i / nframes
    
    n = int(n_slider.val)
    l = int(l_slider.val)
    m = int(m_slider.val)
    
    R_vals = R_nl(R, n, l)
    Y_vals = sph_harm(m, l, THETA, PHI)
    psi = R_vals * Y_vals
    psi_t = psi * np.exp(-1j * t)
    prob_density = np.abs(psi_t)**2
    prob_density /= np.max(prob_density)
    
    X = prob_density * np.sin(THETA) * np.cos(PHI)
    Y = prob_density * np.sin(THETA) * np.sin(PHI)
    Z = prob_density * np.cos(THETA)
    
    surf.remove()
    surf = ax.plot_surface(X, Y, Z, rstride=4, cstride=4, cmap='coolwarm')
    
nframes = 36


anim = FuncAnimation(fig, animate, frames=nframes + 1, interval=2000 / (nframes + 1))
plt.show()

Obviously when executing the code I get no animation since I take the absolute value of the exponential function which just equals one. My result looks like this, which are quite different from the expected results.I think, that my problem is somewhere in calculating the radial part.

I mainly tried different scipy functions for the Laguerre polynomials and spherical harmonics but it didn’t change anything.