This is one way I know of to solve a system of differential equations.
T = ode_solver()
f = lambda t, y: [y[2], y[3], -y[1] - 1, -y[0]]
T.function = f
T.ode_solve(y_0=[1, 1, 0, 0], t_span=[0, 20], num_points=1000)
Just one thing: if the -y[1]-1 and y[0] are replaced with much larger functions of y[0], y[1]. y[2], and y[3], how can I declare those functions before defining T?
