ivp.py

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

a, y0, g0 = 3, 100, 100
T = np.linspace(0, 30, 30)

gs = []
for t in T:
    # G depends on Y
    def g_diff(t, g):
        # Y could be solved iteratively:
        # def y_diff(t, y):
        #     return -y/a
        # solY = solve_ivp(y_diff, [0, t], [y0])
        # y = solY.y[0][-1]

        # But it is just an exponential decay
        y = y0 * np.exp(-t/a)

        return y - min(50, g)

    # Solving G
    solG = solve_ivp(g_diff, [0, t], [g0])

    gs.append(solG.y[0][-1])

# Plotting g & y
plt.plot(T, gs, label='g')
plt.legend()
plt.show()