Accurately solving system of differential equations

Question

So I am trying to solve two equations simultaneously. The goal is to find values for $\frac{de}{dt}$ and $\frac{d}{dt}$ which are the rates of change of the variables $a$ and $e$. I am then progressing the system through time by multiplying the rates by a certain time step and adding this to the starting values to obtain new values for $a$ and $e$.

The system is then solved again and again with these new values for over 1,000,000 data points. The issue is, I know this is not the most accurate method to do this however I am not sure what is. Any help is appreciated! The two equations are shown below:

$$\frac{1}{e}\frac{de}{dt}=-\left[2.0\times10^{21}+ \frac{1.8\times10^{20}}{8} \right] a^{-\frac{13}{2}}$$ $$\frac{1}{a}\frac{da}{dt}=-\left[8.1\times10^{21}e^2+ \frac{2.1\times10^{20}}{8} \right] a^{-\frac{13}{2}}$$

Answer 1

Congratulations, you have rediscovered the Euler method. It is indeed not the most accurate method for solving ordinary differential equations. Moreover, the errors it makes tend to be consistent in direction. Without having tested it for your case, your best choice is probably an adaptive Runge–Kutta method. Such methods should be available in any programming language generally suited for such tasks. Do not try to implement those yourself, unless you really know what you are doing and have a good reason.

As a sidenote, when using such a module, you should normalise your system such that the scale of your dynamical variables ($a$ and $e$) and expected dynamics is roughly $1$ (cf. this answer of mine).