import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def sfr_reheat(y, t, a11, a12, a21, a22, b1, b2, Pd):
omega, z = y
dydt = [a11*omega + a12*z + b1*Pd, a21*omega + a22*z + b2*Pd]
return dydt
H = 5
D = 1
Tr = 8
Fh = 0.3
R = 0.05
Pd = 0.1
a11 = (-Fh/R-D)/(2*H)
a12 = (1 – Fh)/(2*H)
a21 = -1/(R*Tr)
a22 = -1/Tr
b1 = -1/(2*H)
b2 = 0
y0 = [0, 0]
t = np.linspace(0, 30, 300)
sol = odeint(sfr_reheat, y0, t, args=(a11, a12, a21, a22, b1, b2, Pd))
f = sol[:, 0]*50 + 50
plt.plot(t, f)
plt.grid()
plt.show()