Solution (source code)

= Solution

Diagonalize the real <Laplace-Lagrange secular matrix> $\mathbf A$. Its <eigenvalue>s are
$$
\boxed{g_{1,2}=\frac{A_{11}+A_{22}}2
\pm\frac12\sqrt{(A_{11}-A_{22})^2+4A_{12}A_{21}}},
$$
and choose corresponding real <eigenvector>s $\mathbf v_k$. If $V=(\mathbf v_1\ \mathbf v_2)$, the initial <complex eccentricity> vector determines complex mode coefficients
$$
\mathbf c=V^{-1}\mathbf z(0),
\qquad c_k=|c_k|e^{i\beta_k}.
$$
The <matrix exponential> solution is
$$
\mathbf z(t)=\sum_{k=1}^2c_k\mathbf v_k e^{ig_kt}.
$$
Thus
$$
\boxed{z_j(t)=\sum_{k=1}^2e_{jk}e^{i(g_kt+\beta_k)}},
\qquad e_{jk}=|c_k|v_{jk},
$$
where a negative eigenvector component may equivalently be made positive by adding $\pi$ to its phase. The $g_k$ are secular precession frequencies, each eigenvector fixes the planets' eccentricity ratio and relative apsidal orientation, and $\beta_k$ fixes the phase selected by the initial conditions.