Attracteur de Lorenz - RK4

Bonjour
Dans le cadre d'un projet d'étude (3ème année de Licence), je voudrais mettre en place, sur Python, la méthode de RK4 pour résoudre les équations de Lorentz numériquement.
J'ai réussi à le faire sans trop de difficulté avec les méthodes de Newton, mais pour RK4, je bloque !!
Voici ou j'en suis pour le moment.

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]

## Euler NS-Lorenz
def lorenz(x, y, z, s=10, r=28, b=2.667):
    """
    Entrée:
       x, y, z: point dans R3 auquel nous nous interessons
       s, r, b: paramètres de l'attracteur de Lorenz
    Sortie:
       dérivées partiels aux points x, y, z liées aux équations de Lorenz
    """
    x_dot = s*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z
    return x_dot, y_dot, z_dot

pas = 0.01
num_iter = 100

# +1 pour prendre les CI
x_sol = np.empty(num_iter + 1)
y_sol = np.empty(num_iter + 1)
z_sol = np.empty(num_iter + 1)

# CI : Conditions initiales
x_sol[0], y_sol[0], z_sol[0] = (0., 1., 1.05)

# Euler
for i in range(num_iter):
    x_dot, y_dot, z_dot = lorenz(x_sol[ i], y_sol[ i], z_sol[ i])
    x_sol[ i + 1] = x_sol[ i] + (x_dot * pas)
    y_sol[ i + 1] = y_sol[ i] + (y_dot * pas)
    z_sol[ i + 1] = z_sol [ i] + (z_dot * pas)


# Plot
ax = plt.figure().add_subplot(projection='3d')

ax.plot(x_sol, y_sol, z_sol, lw=0.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Attracteur de Lorenz")

plt.show()


## RK4 NS-Lorenz
s,r,b=10,28,2.667
#Lorenz
def f(t,x,y,z):
    return s*(y-x)
def g(t,x,y,z):
    return r*x-y-x*z
def h(t,x,y,z):
    return x*y-b*z

# Initialisation
x,y,z=[0],[1],[1.05]
pas=0.01
num_iter=100 #t=pas*num_iter
t=np.linspace(0,pas*num_iter,num_iter+1)
# RK4

for i in range (num_iter):
    k1=pas*f(t[ i],x[ i],y[ i],z[ i])
    l1=pas*g(t[ i],x[ i],y[ i],z[ i])
    m1=pas*f(t[ i],x[ i],y[ i],z[ i])
    
    k2=f(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))
    l2=g(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))
    m2=h(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))

    k3=f(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))
    l3=g(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))
    m3=h(t[ i]+pas/2,(x[ i]+k1*pas/2),(y[ i]+(l1*pas/2)),(z[ i]+(m1*pas/2)))
    
    k4=f(t[ i]+pas,(x[ i]+k3*pas),(y[ i]+l3*pas),(z[ i]+m3*pas))
    l4=g(t[ i]+pas,(x[ i]+k3*pas),(y[ i]+l3*pas),(z[ i]+m3*pas))
    m4=h(t[ i]+pas,(x[ i]+k3*pas),(y[ i]+l3*pas),(z[ i]+m3*pas))
    
    x.append(x[ i]+pas*(k1+2*k2+2*k3+k4)/6)
    y.append(y[ i]+pas*(l1+2*l2+2*l3+l4)/6)
    z.append(z[ i]+pas*(k1+2*k2+2*k3+k4)/6)

# Plot
ax = plt.figure().add_subplot(projection='3d')

ax.plot(np.array(x), np.array(y), np.array(z), lw=0.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Attracteur de Lorenz")

plt.show()

La première est la figure que j'obtiens avec RK4 (convergence vers l'origine), la seconde avec la méthode de Newton (aspect chaotique). (avec les mêmes constantes et les mêmes conditions initiales). Par conséquent j'imagine qu'il y a une erreur dans mon script d'autant plus que ma version avec RK4 ne varie que très peu si ce n'est pas du tout en fonction des constantes liées au systèmes et des CI.

Merci par avance pour vos retours.
Bien cordialement.
Chaka