Solution (source code)

= Solution

Differentiating the <disturbing function> in <Lagrange planetary equations> gives
$$
\dot e=\sum_j A_j e_j\sin(\varpi-\varpi_j),\qquad
\dot\varpi=A+\sum_j A_j\frac{e_j}{e}\cos(\varpi-\varpi_j).
$$
For the <complex eccentricity> $z=e e^{i\varpi}$,
$$
\dot z=e^{i\varpi}(\dot e+ie\dot\varpi).
$$
Use $\sin\delta+i\cos\delta=i e^{-i\delta}$. Each planetary term becomes $iA_je_je^{i\varpi_j}=iA_jz_j$, so
$$
\boxed{\dot z=iAz+i\sum_{j=1}^N A_jz_j.}
$$
Unlike the polar <orbital element> equations, this linear equation is nonsingular at $e=0$.