Solution (source code)

= Solution

For this <chemical reaction network>, the <law of mass action> gives the reaction rates are $k_1s^2e$, $k_2c$ and $k_3c$. Their <stoichiometry> gives
$$
\dot s=-2k_1s^2e+2k_2c,\quad
\dot e=-k_1s^2e+(k_2+k_3)c,\quad
\dot c=k_1s^2e-(k_2+k_3)c,\quad
\dot p=2k_3c.
$$
Here the convention for $k_1$ absorbs any combinatorial factor for two identical reactants. Direct differentiation gives the two <conservation laws>
$$
\boxed{e+c=e_0,\qquad s+2c+p=s_0.}
$$
Eliminating $e$ and $p$ leaves
$$
\dot s=-2k_1s^2(e_0-c)+2k_2c,\qquad
\dot c=k_1s^2(e_0-c)-(k_2+k_3)c.
$$
For positive $s_0,e_0,k_1$, use the <dimensionless variables> $u=s/s_0$, $v=c/e_0$ and $\tau=k_1e_0s_0t$. With prime marks denoting differentiation in $\tau$, this becomes
$$
\boxed{u'=-2u^2(1-v)+2(\mu-\lambda)v,\qquad\epsilon v'=u^2(1-v)-\mu v,}
$$
where
$$
\boxed{\epsilon=\frac{e_0}{s_0},\quad\mu=\frac{k_2+k_3}{k_1s_0^2},\quad\lambda=\frac{k_3}{k_1s_0^2},\quad u(0)=1,\quad v(0)=0.}
$$