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}')