Applying Forward Euler Method to a Three-Box Model System of ODEs

I am trying to model a system of coupled ODEs which represent a three-box ocean model of phosphorous concentration (y) in the low-latitude surface ocean (Box 1), high-latitude deep ocean (Box 2), and deep ocean (Box 3). The ODEs are given below:

dy1dt = (Q/V1)*(y3-y1) - y1/tau                                        # dP/dt in Box 1
dy2dt = (Q/V2)*(y1-y2) + (qh/V2)*(y3-y2) - (fh/V2)                     # dP/dt in Box 2
dy3dt = (Q/V3)*(y2-y3) + (qh/V3)*(y2-y3) + (V1*y1)/(V3*tau) + (fh/V3)  # dP/dt in Box 3

The constants and box volumes are given by:

### Define Constants
tau = 86400              # s
VT  = 1.37e18            # m3 
Q   = 25e6               # m3/s 
qh  = 38e6               # m3/s
fh  = 0.0022             # mol/m3
avp = 0.00215            # mol/m3
    
### Calculate Surface Area of Ocean
r = 6.4e6                # m
earth = 4*np.pi*(r**2)   # m2
ocean = .70*earth        # m2   
    
### Calculate Volumes for Each Box
V1 = .85*100*ocean       # m3
V2 = .15*250*ocean       # m3
V3 = VT-V1-V2            # m3

This can be put into matrix form y = Ay + f, where y = [y1, y2, y3]. I have provided the matrices and initial conditions below:

A = np.array([[(-Q/V1)-(1/tau),0,(Q/V1)],
                  [(Q/V2),(-Q/V2)-(qh/V2),(qh/V2)],
                  [(V1/(V3*tau)),((Q+qh)/V3),((-Q-qh)/V3)]])
f = np.array([[0],[-fh/V2],[fh/V3]])

y1 = y2 = y3 = 0.00215 # mol/m3

I am having trouble adapting the Forward Euler method to apply to a system on linear ODEs, rather than just one. This is what I have come up with so far (it runs with no issues but doesn't work if that makes sense; I think it has something t do with initial conditions?):

### Define a Function for the Euler-Forward Scheme
import numpy as np
def ForwardEuler(t0,y0,N,dt):
        N = 100
        dt = 0.1
        # Create empty 2D arrays for t and y
        t = np.zeros([N+1,3,3]) # steps, # variables, # solutions
        y = np.zeros([N+1,3,3])
        # Assign each ODE to a vector element
        y1 = y[0]
        y2 = y[1]
        y3 = y[2]
        # Set initial conditions for each solution
        t[0, 0] = t0[0] 
        y[0, 0] = y0[0]
        t[0, 1] = t0[1] 
        y[0, 1] = y0[1]
        t[0, 2] = t0[2] 
        y[0, 2] = y0[2]
        
        for i in trange(int(N)):
            t[i+1]  = t[i]  + t[i]*dt
            y1[i+1] = y1[i] + ((Q/V1)*(y3[i]-y1[i]) - (y1[i]/tau))*dt
            y2[i+1] = y2[i] + ((Q/V2)*(y1[i]-y2[i]) + (qh/V2)*(y3[i]-y2[i]) - (fh/V2))*dt
            y3[i+1] = y3[i] + ((Q/V3)*(y2[i]-y3[i]) + (qh/V3)*(y2[i]-y3[i]) + (V1*y1[i])/(V3*tau) + (fh/V3))*dt
        return t, y1, y2, y3

Any help on this is greatly appreciated. I have not found any resources online that go through the Euler Forward for a system of 3 ODEs, and am at a loss. I am happy to explain further if there are more questions.

Solution

As Lutz Lehmann pointed out, you need to design a simple ODE system. You could define the whole ODE system inside a function as follows:

import numpy as np

def fun(t, RHS):

    # get initial boundary condition values
    y1 = RHS[0]
    y2 = RHS[1]
    y3 = RHS[2]

    # calculte rate of respective variables
    dy1dt = (Q/V1)*(y3-y1) - y1/tau
    dy2dt = (Q/V2)*(y1-y2) + (qh/V2)*(y3-y2) - (fh/V2)
    dy3dt = (Q/V3)*(y2-y3) + (qh/V3)*(y2-y3) + (V1*y1)/(V3*tau) + (fh/V3)

    # Left-hand side of ODE
    LHS = np.zeros([3,])

    LHS[0] = dy1dt
    LHS[1] = dy2dt
    LHS[2] = dy3dt

    return LHS

In the above function, we get time t as an argument and the initial values of y1, y2, and y3 as a list in the variable RHS, which is then unpacked to get the respective variables. Afterward, the rate equation of each variable is defined. In the end, the calculated rates are returned also as a list in the variable LHS.

Now, we can define a simple Euler Forward method to solve this ODE system as follows:

def ForwardEuler(fun,t0,y0,N,dt):

    t = np.arange(t0,N+dt,dt)
    y = np.zeros([len(t), len(y0)])
    y[0] = y0

    for i in range(len(t)-1):
        y[i+1,:] = y[i,:] + fun(t[i],y[i,:]) * dt

    return t, y

Here, we create a time range from 0 to 100 with a step size of 0.1. Afterward, we create an array of zeros with the shape (len(t), len(y0)) which is in this case (1001,3). We need to do this because we want to solve fun for the range of t (1001) and the RHS variable of fun has a shape of (3,) ([y1, y2, y3]). So for each and every point in t, we will solve for the three variables of RHS, which will be returned as LHS.

In the end, we can solve this ODE system as follows:

dt = 0.1
N = 100
y0 = [0.00215, 0.00215, 0.00215]
t0 = 0

t,y = ForwardEuler(fun,t0,y0,N,dt)

Solution using scipy.integrate

As Lutz Lehmann also pointed out, you can use scipy.integrate for this purpose as well which is far easier. Here you can use the above defined fun and simply solve the ODE as follows:

import numpy as np
from scipy.integrate import odeint

dt = 0.1
N = 100
t = np.linspace(0,N,int(N/dt))
y0 = [0.00215, 0.00215, 0.00215]

res = odeint(fun, y0, t, tfirst=True)
print(res)