A mini solar system

The presence of other planets implies that the total force on a planet is no longer a central force. Furthermore, since the orbits are not exactly on the same plane, the analysis must be extended to 3D. However, for simplicity, we are going to consider a two-dimensional solar system, with two planets in orbit around the sun.

The equations of motion of the two planets of mass $m_1$ and $m_2$ can be written in vector form as

$$m_1\frac{d^2 {\mathbf r}_1}{dt^2}=-\frac{m_1MG}{r_1^3}{\mathbf r}_1+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21},$$ (33)

$$m_2\frac{d^2 {\mathbf r}_2}{dt^2}=-\frac{m_2MG}{r_2^3}{\mathbf r}_2+\frac{m_1m_2G}{r_{21}^3}{\mathbf r}_{21},$$ (34)

where ${\mathbf r}_1$ and ${\mathbf r}_2$ are directed form the sun to the planets, and ${\mathbf r}_{21}={\mathbf r}_2-{\mathbf r}_1$ is the vector from planet 1 to planet 2. This is a problem with no analytical solution, but its numerical solution can be obtained extending our previous analysis for the two-body problem.