Let the two states be $\lambda:[0,T]\rightarrow \mathbb{R}^+$ and $p:[0,T]\rightarrow\mathbb{R}^{++}$. $b$ is a constant which is positive.
We have $$ \begin{align} \dot{\lambda} &=-\frac{b}{\lambda}\log (pb) \\ \dot{p} &=[\frac{1}{\lambda^2}-\frac{\log (pb)}{\lambda^2}]pb \end{align} $$ with $\lambda(0)=\lambda_0>0$ and $p(T)=p_T>0$.
How to numerically solve this nonlinear differential equations with the boundary conditions mentioned above? Sorry if there is any misuse of terms. It would be great if you guys can provide some materials about how to solve this kind of problem.