Solution (source code)

= Solution

For the relative position $\mathbf r=\mathbf r_1-\mathbf r_2$, <Newton's law of universal gravitation> gives the <two-body problem> in its reduced one-body form,
$$
\ddot{\mathbf r}=-\frac{\mu}{r^3}\mathbf r,
\qquad
\mu=G(M_1+M_2).
$$
Taking the <dot product> with the relative velocity $\mathbf v=\dot{\mathbf r}$ gives
$$
\frac d{dt}\left(\frac12v^2\right)
=-\frac{\mu}{r^3}\mathbf r\mathbin{\cdot}\dot{\mathbf r}
=\frac d{dt}\left(\frac\mu r\right).
$$
Integration therefore yields conservation of <specific orbital energy>:
$$
\boxed{\frac12v^2-\frac\mu r=C}.
$$