I've been trying to implement a higher order Runge-Kutta numerical integrator in [MATLAB](https://www.mathworks.com/products/matlab.html) for a bit now with ambitions of n-body trajectory fun. Unfortunately, I have not had success in reaping the rewards of such higher order methods.

I have (I believe) a good grasp of how Runge-Kutta methods work and have ample experience with lower order methods like an [RK4](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method) or [Dormand Prince](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method).

Various answers on the site have suggested [RK78](https://space.stackexchange.com/a/55922/40257) and [RK89](https://space.stackexchange.com/a/55832/40257). I found [Fehlberg's Technical Report (NTRS ID: 19680027281)](https://ntrs.nasa.gov/citations/19680027281) describing both, as well as C++ source code from [Trick](https://github.com/nasa/trick/blob/master/trick_source/er7_utils/integration/rkf78/src/rkf78_butcher_tableau.cc) and [GMAT](https://github.com/ChristopherRabotin/GMAT/blob/GMAT-2020a/src/base/propagator/RungeKutta89.cpp) giving the [Butcher tableaus](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#:%7E:text=These%20data%20are%20usually%20arranged%20in%20a%20mnemonic%20device%2C%20known%20as%20a%20Butcher%20tableau) (in a C++ structure) for each method. This gives me confidence that my coefficients are correct as I copied them into my MATLAB script (and only "translated" to MATLAB).

My simple test case is an object in a 2D 185 km orbit above "Earth", that is, a point mass with no atmosphere. I used a fixed 30 second time step (just trying to get *one* of the orders to work!).

For both of my RK78 & RK89 implementations I get nearly identical erroneous outputs of this growing spiral-like trajectory:

---

**Edit:** my specific energy plots were not correct, I have fixed those, the problem(s) still remain with the trajectory.

---

RK78:

[![RK78 results](https://i.stack.imgur.com/cn9vo.png)](https://i.stack.imgur.com/cn9vo.png)

RK89:

[![RK89 results](https://i.stack.imgur.com/hN3ml.png)](https://i.stack.imgur.com/hN3ml.png)

These results occur regardless of whether I select the 7th or 8th (or 8th or 9th) order output.

Compare these to a much more reasonable result from a [Dormand Price](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method) integration programmed in the same structure (trajectory plot excluded because it's expectedly benign):

[![RK45 specific energy](https://i.stack.imgur.com/11mZJm.png)](https://i.stack.imgur.com/11mZJm.png)

Here is a code snippet of the core iteration scheme:

```
% k is a m x 4 matrix where m is the number of stages, 16 for RK89 (shown), 13 for RK78.
% multiplying the weights (C, m x 1, transposed) by k (m x 4) gives the output 'slopes' (1 x 4)

while 1
    % Runge-Kutta Slopes
    for j = 1:16 % k_j (individual k slopes)
        input = zi(i,:);
        for jj = 1:16 % (component of k slopes, B matrix)
           input = input + B(j,jj)*k(jj,:);
        end
        k(j,:) = dz(input,G,M_E)'; % call physics model for j'th node
    end

    % Forward local extrapolation
    i = i+1; % move time forward

    % Break Conditions
    ...

    % Integrate
    zi(i,:) = zi(i-1,:) + dt*(C'*k); % next state
end

% force model function 'dz':
function dzdt = dz(z,G,M_E)
aG_E = G*M_E/(z(1)^2 + z(2)^2); % gravitational acceleration (m/s^2)

% State Variable Definitions:
% z1 = x, z2 = y,
% z3 = Vx, z4 = Vy

% 'slopes':
dzdt(1) = z(3);
dzdt(2) = z(4);
dzdt(3) = -aG_E*z(1)/sqrt(z(1)^2 + z(2)^2);
dzdt(4) = -aG_E*z(2)/sqrt(z(1)^2 + z(2)^2);
end

```

**Does anyone have any clues as to what is going (*so horribly*) wrong?**