<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;"># xâ€˜= Ïƒ(y â€“ x)
# yâ€˜ = x(Ï â€“ z) â€“ y
# zâ€˜ = xy â€“ Î²z


import numpy
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def lorenz(t, xyz):
    x, y, z = xyz
    s, r, b = 10, 28, 8 / 3.0  # parameters Lorentz used
    return [s * (y - x), x * (r - z) - y, x * y - b * z]


a, b = 0, 40
t = numpy.linspace(a, b, 4000)

sol1 = solve_ivp(lorenz, [a, b], [1, 1, 1], t_eval=t)
sol2 = solve_ivp(lorenz, [a, b], [1, 1, 1.00001], t_eval=t)

plt.plot(sol1.y[0], sol1.y[2])
plt.xlabel("$x$")
plt.ylabel("$z$")
plt.savefig("lorenz_xz.png")
plt.show()
plt.close()


plt.subplot(211)
plt.plot(sol1.t, sol1.y[0])
plt.xlabel("$t$")
plt.ylabel("$x_1(t)$")
plt.subplot(212)
plt.plot(sol1.t, sol1.y[0] - sol2.y[0])
plt.ylabel("$x_1(t) - x_2(t)$")
plt.xlabel("$t$")
plt.savefig("lorenz_x.png")
plt.show()
</pre></body></html>