Solution (source code)

= Solution

Multiply the <master equation> by $z^n$ and sum over $n$. Shifting indices produces
$$
\boxed{\partial_tF=(z-1)\sum_nq_nz^nP_n+(z^{-1}-1)\sum_nr_nz^nP_n.}
$$
The apparent negative power causes no problem at $z=0$ when $r_0=0$. Let $D=z\partial_z$, so $Dz^n=nz^n$. If the rates are polynomials in $n$, replace their arguments by $D$ to obtain a finite-order <partial differential equation>:
$$
\boxed{\partial_tF=\bigl[(z-1)q(D)+(z^{-1}-1)r(D)\bigr]F.}
$$
For arbitrary rates the weighted sums still give an exact generating-function evolution, but need not reduce to a finite-order differential operator. The <probability generating function> obeys
$$
\boxed{F(1,t)=1,}
$$
which is the normalization condition. The evolution has $\partial_tF(1,t)=0$ when the rate sums converge, so an initially normalized distribution stays normalized. Its derivatives at $z=1$ also give the factorial <moments>, and expansion around $z=0$ recovers individual <probabilities>.