Here is a Lorenz ODE solver in 20 lines of dependency-free Python to enlighten you! Find out more about the Taylor Series Method at https://github.com/m4r35n357/ODE-Playground, read the docs and follow the references, marvel at how long this technique has lain undiscovered by the masses.
#!/usr/bin/env python3
from sys import argv
# Example: ./tsm.py lorenz 10 .01 3000 -15.8 -17.48 35.64 10 28 8 3
model, order, δt, n_steps = argv[1], int(argv[2]), float(argv[3]), int(argv[4]) # integrator controls
x0, y0, z0 = float(argv[5]), float(argv[6]), float(argv[7]) # initial values
x, y, z = [0.0] * (order + 1), [0.0] * (order + 1), [0.0] * (order + 1) # Taylor coefficient jets
σ, ρ, β = float(argv[8]), float(argv[9]), float(argv[10]) / float(argv[11]) # Lorenz parameters
print(f'{x0:.9e} {y0:.9e} {z0:.9e} {0.0:.5e}')
for step in range(1, n_steps + 1):
x[0], y[0], z[0] = x0, y0, z0 # build the coefficient jets
for k in range(order):
x[k + 1] = σ * (y[k] - x[k]) / (k + 1)
y[k + 1] = (ρ * x[k] - sum(x[j] * z[k - j] for j in range(k + 1)) - y[k]) / (k + 1)
z[k + 1] = (sum(x[j] * y[k - j] for j in range(k + 1)) - β * z[k]) / (k + 1)
x0, y0, z0 = x[-1], y[-1], z[-1] # perform an integration step using Horner's method
for i in range(len(x) - 2, -1, -1):
x0 = x0 * δt + x[i]
y0 = y0 * δt + y[i]
z0 = z0 * δt + z[i]
print(f'{x0:.9e} {y0:.9e} {z0:.9e} {step * δt:.5e}')