Angular position as a function of time for elliptical orbits 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.
 A: 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.
A: this is the differential equation that you want to solve
$$\dot\theta=2\,{\frac {\pi \, \left( 1+e\cos \left( \theta \right)  \right) ^{2}}{
P \left( 1-{e}^{2} \right) ^{3/2}}}
=f(\theta)$$
you can use the explicit  Euler method
$$\theta_{t+h}=\theta_t+\left[1-\frac{\partial F}{\partial \theta}\bigg|_{\theta_t}\right]^{-1}\,F(\theta_t)$$
where $~F=h\,f(\theta)~$ , h is the step time
and t is the time. the advantech of this method is that you can obtain the derivative   $~\frac{\partial F}{\partial \theta}~$ analytical
result

