from scipy.integrate import solve_ivp
import numpy as np
k1, k2 = 1, 1
def reaction_model(t, y):
a, b, c = y
dadt = -k1 * a
dbdt = k1 * a - k2 * b
dcdt = k2 * b
return [dadt, dbdt, dcdt]
t_span = 10
t_eval = np.linspace(0, t_span, 100)
iv = [1, 0, 0]
sol = solve_ivp(reaction_model, t_span=[0, t_span], t_eval=t_eval, y0=iv)
