The following system of coupled ODEs arises in the study of DNA sequence evolution:

\begin{eqnarray*}
  \frac{da}{dt} & = & \frac{\mu (1-y) b u}{S - y(S-b-v)} - (\lambda +\mu ) a \\
  \frac{db}{dt} & = & -\frac{\mu (b+v) b}{S - y(S-b-v)} + \lambda (1-b) \\
  \frac{du}{dt} & = & -\frac{\mu (b+v) u}{S - y(S-b-v)} + \lambda a \\
  \frac{dv}{dt} & = & \frac{\mu (b+v) (S-v)}{S - y(S-b-v)} \\
  S & = & \exp\left(\frac{\lambda t}{1-x}\right)-1
\end{eqnarray*}
for $t>0$, with $a(0)=1$, $b(0)=u(0)=v(0)=0$,
$a'(0)=-\lambda-\mu$, $b'(0)=u'(0)=\lambda$, and $v'(0)=0$.

Specifically, they are a moment-matched approximation to the General Geometric Indel model, the simplest continuous-time Markov chain on strings that allows insertions and deletions (a.k.a. "indels") of length > 1; details are [here](https://academic.oup.com/genetics/article/216/4/1187/6065876).

The parameters are $\lambda,\mu \in \Re^+$, the indel rates, and $x,y \in [0,1)$, the corresponding probability parameters for the geometric distributions over indel lengths. The mean insertion length is $1/(1-x)$ and the mean deletion length is $1/(1-y)$.

The special case $x=y=0$ (where indels all have unit length) is known as the [TKF91 model](https://www.southgreen.fr/glossary/term/682), and has the following exact solution:

\begin{eqnarray*}
  a & = & \frac{\mu - \lambda}{\mu L - \lambda M} LM \\
  b & = & \frac{\lambda}{\mu L - \lambda M} (L - M) \\
  u & = & \frac{\mu - \lambda}{\mu L - \lambda M} M(1-L) \\
  v & = & \frac{\mu - \lambda}{\mu L - \lambda M} - 1
\end{eqnarray*}
with $L=\exp(-\lambda t)$ and $M=\exp(-\mu t)$.

Can anyone suggest a way to solve this? Things I have tried/observed:

 1. Clearly $(b,v)$ form an independent subsystem of coupled ODEs, so it's presumably worth solving these first.
 2. I haven't yet been able to simplify further to just a single ODE.
 3. I've tried plugging into Mathematica's DSolve as-is. No joy.
 4. Following the TKF91 solutions, I tried using a trial solution that is a rational function of polynomials that are first-order in $L=\exp(-\lambda t/(1-x))$ and $M=\exp(-\mu t/(1-y))$ and share a common denominator, i.e. $b=\frac{b_0 + b_L L + b_M M + b_{LM} LM}{d_0 + d_L L + d_M M + d_{LM} LM}$ and $v=\frac{v_0 + v_L L + v_M M + v_{LM} LM}{d_0 + d_L L + d_M M + d_{LM} LM}$ (Mathematica gets stuck; if I constrain $d_0=d_{LM}=0$, as in TKF91, it says there is no solution.)
  5. I can of course approximate it numerically, e.g. with Runge Kutta methods, but I'd really like an analytical solution.