odeint.py

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint

def Lorenz(X,t):
    x, y, z = X    
    dx = -σ * x + σ * y
    dy = ρ * x - y - x * z
    dz = -β * z + x * y
    return [dx, dy, dz]
    
x0 = 0.0
y0 = 1.0
z0 = 1.05

σ  = 10.0
ρ  = 28.0
β  = 8.0/3.0

tlm  = [0.0,30.0]

H  = [1e-3,1e-4,1e-5,1e-6]

T = []
Y = []
for i,h in enumerate(H):
    print(i)
    
    t = np.arange(tlm[0],tlm[1],h)

    y = odeint(Lorenz,[x0,y0,z0],t)

    Y.append(y)
    T.append(t)

plt.close('all')
fig = plt.figure(figsize=(6,6))

for i, y in enumerate(Y):
    ax = fig.add_subplot(2, 2, i+1, projection='3d')
    ax.plot(y[:,0],y[:,1],y[:,2],lw=0.5)
    plt.title("%.1e" % H[i])
plt.savefig('Lorenz_odeint.png')

plt.figure()
for i, y in enumerate(Y):
    plt.plot(T[i],y[:,0],label="%.1e" % H[i],linewidth=1)
plt.legend(fontsize=8,loc=1)
plt.xlabel('t')
plt.ylabel('x')
plt.grid()
plt.tight_layout()
plt.savefig('Lorenz_odeint2.png')

plt.figure()
for i, y in enumerate(Y):
    plt.plot(T[i],y[:,0],label="%.1e" % H[i],linewidth=1)
plt.legend(fontsize=8,loc=1)
plt.grid()
plt.xlabel('t')
plt.ylabel('x')
plt.xlim([20,30])
plt.tight_layout()
plt.savefig('Lorenz_odeint3.png')