The differential equation (REF ) governing \(\phi _1\)  is more difficult due to the presence of a mutation term involving \(\phi _2\) .
We now have a single inhomogeneous nonlinear Riccati equation that cannot be solved by a simple substitution trick.
However, we note that it is straightforward to solve \(\phi _1\)  numerically using any number of differential equation solvers; we use DifferentialEquations.jl package in the Julia programming language [1]}, [2]}. The computational bottleneck to evaluate the joint PGF at any time \(t>0\)  amounts to numerically solving a single ODE.
<FIGURE>