It is a well-known fact that there is no analytical solution to the two-body problem as a function of time. However, I wanted to derive the differential equation for the angular position of an elliptical orbit, such that I could use it to code a numerical solution.
According to Kepler's laws,
$$\frac{1}{2}r^2\frac{\mathrm{d}\theta}{\mathrm{d}t} = \frac{\pi a b}{P}$$
Where $a$ and $b$ are the semi-major and semi-minor axis.
Since
$$r=\frac{a(1-e^2)}{1+e \cos(\theta)}$$
and
$$b=a \sqrt{1-e^2}$$
it means that
$$\frac{\mathrm{d}\theta}{\mathrm{d}t} = \frac{2\pi a^2 \sqrt{1-e^2}}{Pr^2}$$
$$\frac{\mathrm{d}\theta}{\mathrm{d}t} = \frac{2\pi(1+e \cos(\theta))^2 }{P (1-e^2)^{3/2}} $$
Which, according to WolframAlpha, has a solution for $t$ in terms of $\theta$. In this sense, I wanted to know if this equation is correct and, if so, suggestions of methods for computing a numerical solution. I have tried the well-known Euler's method, but I don't think it is the most proper way.
Best Answer
Your final differential equation for $\theta(t)$ looks correct.
Euler's method for solving a first-order differential equation is the easiest method. There are other methods (see Numerical methods for ordinary differential equations - Methods) which are more elaborated and give more precise results. I suggest to try a Runge-Kutta method.
For nearly circular orbits (i.e. for $e\ll 1$) there is also the option to use a Fourier expansion for $\theta(t)$. Taken from Mean anomaly - Formula: $$\begin{align} \theta(t) &= M \\ &+ \left(2e - \frac{1}{4}e^3 + \cdots \right)\sin(M) \\ &+ \left(\frac{5}{4}e^2 + \cdots \right) \sin(2M)\\ &+ \left(\frac{13}{12}e ^3 + \cdots \right) \sin(3M) \\ &+ \cdots \end{align}$$ where $M=\frac{2\pi}{P}t$ is the so-called mean anomaly. And $\cdots$ stands for higher-order terms in $e$ which become neglectable in case of nearly circular orbits.