How to change this Runge-Kutta method implementation from first order ode solver to system of ODEs solver?

Question

I implemented 4-step Runge-Kutta method (k1..k4) ODE solver for a function $u'(x) = f(x,u(x))$ with initial condition $u(x_0) = u_0$

But it solves just ODEs of the first order. How could I change the code, so that it also takes ODEs of higher order? ${\bf u}'(x)={\bf f}(x,{\bf u}(x))$ with initial condition ${\bf u}(x_0)= {\bf u}_0$? All as vectors here not scalars.

Here is the code for my solution, which takes first order ODEs:

def euler(f, x, y, h):
    yn = y + h*f(x,y)
    xn = x + h
    return xn, yn

def rk4(f, x, y, h):
    k1 = h*f(x,y)
    k2 = h*f(x+(1/2)*h,y+(1/2)*k1) 
    k3 = h*f(x+(1/2)*h,y+(1/2)*k2)
    k4 = h*f(x+h,y+k3)
    yn = y + ((1/6)*(k1+(2*k2)+(2*k3)+k4))
    xn = x + h
    return xn, yn

def integrate(method, f, t0, y0, tend, h): # Depending on the 'method' (euler or rk4) this method solves the ode f)
    T = [t0]
    Y = [y0]
    t = t0
    y = y0 
    while (t < tend):
        h = min(h, tend-t)
        t, y = method(f, t, y, h)
        T.append(t)
        Y.append(y)
    return np.array(T), np.array(Y)

#def f1(t, y):   # this is an 1st order ODE
    #return (t*y) + t

# Usage example
xv, yv = integrate(method=rk4, f=f1, t0=0, y0=2, tend=1, h=0.5)
yv[-1] # Output: 6.15203857421875

Is there a solution, to make this work also for a given system of ODEs? (dynamically, without knowing the input ODE if it's 1st or higher order)for example this one: $$u''(x)+u(x)+u(x)^3 =0 ~~{\rm with}~~ u(0)=u'(0)=0 $$ or this one: $$u_1'(x) = 98u_1(x) + 198u_2(x) ~~{\rm and}~~ u_2'(x) = −99u_1(x) − 199u_2(x) \\~~{\rm with}~~ u_1(0) = 1 ~~{\rm and}~~ u_2(0) = 0$$

EDIT:

changed the code to this now, to try if the 1st given example 2nd order ODE works:

def f(t,y):
    return np.array([ y[1], -y[0]-(y[0])**3 ])

Y0 = np.array([0.0,0.0])

xv, yv = integrate(method=rk4, f=f, t0=0.0, y0=Y0, tend=100.0, h=h)
print("[ %20.15f, %20.15f]"%(yv[-1,0], yv[-1,1]) )
print(yv)

but I get an output consisting of 0 vectors like this:

[    0.000000000000000,    0.000000000000000]
[[0. 0.]
 [0. 0.]
 [0. 0.]
 ...
 [0. 0.]
 [0. 0.]
 [0. 0.]]

Answer 1

There is not much to change, depending on how "user friendly" you want to make it. The quick solution is to explicitly use np.array for vectors, for the classical example of the physical pendulum $y''+\sin(y)=0$ you could code

def f(t, y): return np.array([ y[1], -np.sin(y[0]) ])

y0 = np.array([1.0,0.0])

for h in [0.5, 0.05, 0.005, 0.0005]:
    xv, yv = integrate(method=rk4, f=f, t0=0.0, y0=y0, tend=1.0, h=h)
    print("[ %20.15f, %20.15f]"%(yv[-1,0], yv[-1,1]) )

gives the output

[    0.600232949607472,   -0.754659515900116]
[    0.600085385678949,   -0.754963692069849]
[    0.600085366129492,   -0.754963713951030]
[    0.600085366127465,   -0.754963713953161]

correctly gaining 4 digits for every reduction of the step size by 10.

---

If you wanted to remove the np.array from the user side, you would have to change the integrate function to detect the type of y0 and if it is not a scalar, redefine initial vector and the input function

if y0 is not np.float:
    y0=np.asarray(y0);
    f_=f; f=lambda t,y:np.asarray(f_(t,y));

or similar.