I can solve the following differential equation,
\begin{align} i\frac{d b_1(t)}{dt} &= \alpha b_1(t) + \beta_1 \left[E_g(t)\right] b_2(t)\notag \\ &+ \beta_2 \left[E_p(t)\right] b_2(t)\\ i\frac{d b_2(t)}{dt} &= \alpha_2 b_2(t) + \beta_3 \left[E_g^*(t)\right] b_1(t) \notag \\ &+ \beta_4\left[E_p^*(t)\right] b_1(t) \end{align}
where $E_g(t)$ and $E_p(t)$ are some complex functions of t. This differential equation can be solved numerically to find the functional dependence of $b_1(t)$ and $b_2(t)$. Then, I did, $$P(t) = b_2(t)b_1^*(t)$$ I want to numerically find, \begin{equation} \left. \frac{\partial P(t)}{\partial E_g(t)}\right\vert_{E_g=0,E_p=0} \end{equation} Can someone suggest an algorithm to implement it?