I want to solve a two-body problem using RK2. The following is a system of equations to be solved \begin{equation}\label{two-body} x^{\prime \prime}=-\frac{x}{r^{3}}, y^{\prime \prime}=-\frac{y}{r^{3}}, r=\sqrt{x^{2}+y^{2}} \end{equation}
let \begin{equation} \begin{aligned} v_{1}=x \\ v_{2}=y \\ v_{3}=x^{\prime} \\ v_{4}=y^{\prime} \end{aligned} \end{equation} So we get \begin{equation} \begin{aligned} v_{1}^{\prime}=v_{3} \\ v_{2}^{\prime}=v_{4} \\ v_{3}^{\prime}=-\frac{v_{1}}{r^{3}} \\ v_{4}^{\prime}=-\frac{v_{2}}{r^{3}} \end{aligned} \end{equation} Then let $\textbf{v}=\left[v_1,v_2,v_3,v_4\right]=\left[x, y, x^{\prime}, y^{\prime}\right]$ and $f(t, v)=v^{\prime}(t)=\left[x^{\prime}, y^{\prime},-\frac{x}{r^{3}},-\frac{y}{r^{3}}\right], r=\sqrt{x^{2}+y^{2}}$.
Next, I arrange IVP like this \begin{equation}\label{MSA two body} \begin{aligned} \textbf{v}^{\prime}(t)=f(t, \textbf{v})\\ \textbf{v}\left(t_{0}\right)=\textbf{v}_{0}. \end{aligned} \end{equation}
Now, I want to solve with RK2
$ \begin{aligned} \mathbf{v}_{n+1}=& \mathbf{v}_{n}+h\left(b_{1} k_{1}+b_{2} k_{2}\right) \\ k_{1}=& f\left(t_{n}, \mathbf{v}_{n}\right)=\left[x_{n}^{\prime}, y_{n}^{\prime}, x_{n}^{\prime \prime}, y_{n}^{\prime \prime}\right] \\ k_{2}=& f\left(t_{n}+a_{21} h_{1}, \mathbf{v}_{n}+a_{21} h k_{1}\right) \\ =& {\left[x_{n}^{\prime}+a_{21} h\left(\frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}\right), y_{n}^{\prime}+a_{21} h\left(\frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}\right),\right.} \\ &\left.\frac{-\left(x_{n}+a_{21} h x_{n}^{\prime}\right)}{\left(\left(x_{n}+a_{21} h x_{n}^{\prime}\right)^{2}+\left(y_{n}+a_{21} h y_{n}^{\prime}\right)^{2}\right)^{\frac{3}{2}}}, \frac{-\left(y_{n}+a_{21} h y_{n}^{\prime}\right)}{\left(\left(x_{n}+a_{21} h x_{n}^{\prime}\right)^{2}+\left(y_{n}+a_{21} h y_{n}^{\prime}\right)^{2}\right)^{\frac{3}{2}}}\right] \\ =& {\left[k_{1 x}+a_{21} h k_{1 x^{\prime}}, k_{1 y}+a_{21} h k_{1 y^{\prime}},\right.} \\ &\left.f_{x^{\prime}}\left(x_{n}+a_{21} h k_{1 x}, y+a_{21} h k_{1 y}\right), f_{y^{\prime}}\left(x_{n}+a_{21} h k_{1 x}, y_{n}+a_{21} h k_{1 y}\right)\right] . \end{aligned} $
with $k_{1 x}= x_{n}, k_{1 y}= y_{n}, k_{1 x^{\prime}}= \frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}, k_{1 y^{\prime}}= \frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}$ and $f_{x^{\prime}}\left(x_{n}, y_{n}\right)= \frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}},f_{y^{\prime}}\left(x_{n}, y_{n}\right)= \frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}$.
Now, I will compute analytics and numeric value (with RK2 there are Heun, Ralston,midpoint and new method with unusual RK2 coefficient) from two-body problem then compare the error. Here the program
https://colab.research.google.com/drive/1P0aRbqmFa1b77mGtzsojmZwe3BThzf_D?usp=sharing
And I got plot like this
The question is how can when h decreasing the error gets bigger in all methods after N>600? (N is number of point on interval)