# Accurately solving system of differential equations

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}}$$

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).