 Solution

The variables $y_j=I_j{e}^{i{\Omega }_j}$ are complex inclinations, despite the question's description as eccentricities. Diagonalize the Laplace-Lagrange inclination matrix: let real eigenvectors ${\mathbf{v}}_k$ satisfy $B{\mathbf{v}}_k={\lambda }_k{\mathbf{v}}_k$, and form $V=({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)$. The physical matrix is diagonalizable by positive diagonal symmetrization of a matrix, as shown below.

The matrix exponential solution isWrite $c_k=I_k{e}^{i{\gamma }_k}$ with $I_k\ge 0$, and define $I_{jk}=I_k{v}_{jk}$. This givesThe eigenvalues determine nodal precession frequencies, the eigenvectors determine relative inclinations within each secular eigenmode, and the initial complex inclinations determine their amplitudes and phases. The real $I_{jk}$ can have either sign, encoding opposite nodal phases. The zero-frequency mode is a common fixed tilt.