Example coupled ODEs

[RemDelaporteMathurin]

from scipy.optimize import fsolve

lamb = 1
lamb2 = 0.5
delta = 0.1

c_old = [10, 20]
t = 0
times = []
c1s = []
c2s = []

while t < 1:

    def equation(c):
        c1, c2 = c
        c1_old, c2_old = c_old

        return [
            (c1 - c1_old) / delta - (-lamb * c1),
            (c2 - c2_old) / delta - (-lamb * c2) - lamb2 * c1,
        ]

    new_c = fsolve(equation, x0=[0, 0])

    c_old = new_c
    t += delta
    times.append(t)
    c1, c2 = new_c
    c1s.append(c1)
    c2s.append(c2)

import matplotlib.pyplot as plt

plt.plot(times, c1s)
plt.plot(times, c2s)
plt.show()
``