Solution (source code)

= Solution

Multiplying the tangential <equation of motion> by $r$ gives
$$
\frac d{dt}(r^2\dot\theta)=r^2\ddot\theta+2r\dot r\dot\theta=0.
$$
Therefore $h=r^2\dot\theta$ is constant. It is the signed <specific angular momentum>; reverse the angular orientation if necessary to take $h>0$ for a nonradial <orbit>. Set $u=1/r$ and denote differentiation with respect to $\theta$ by a prime. Then
$$
\dot\theta=hu^2,\qquad \dot r=-hu',\qquad \ddot r=-h^2u^2u'',\qquad r\dot\theta^2=h^2u^3.
$$
Substitution into the radial <equation of motion> and division by $-h^2u^2$ give the <Binet equation>
$$
u''+u=\frac{GM}{h^2}.
$$
Its general solution is $u=GM/h^2+A\cos\theta+B\sin\theta$. Write $A=(GM/h^2)e\cos\theta_0$ and $B=(GM/h^2)e\sin\theta_0$, where $e\geq0$. This produces the <Kepler orbit>
$$
\boxed{\frac{h^2u}{GM}=1+e\cos(\theta-\theta_0).}
$$
Here $e$ is the <orbital eccentricity>, and $\theta_0$ points toward <periapsis> when $e>0$. Only portions with $u>0$ describe positive radial distances. The expression assumes $h\ne0$; a zero-<angular momentum> trajectory is radial and must be treated separately.