To make my comments more explicit:
Try with the following code for the DDE function that eliminates P
function dydt = dde(t,y,yd,A)
n = length(A);
dydt = diag(A).*y;
k = 0;
for i = 1:n
for j = (i+1):n
dydt(i) += A(i,j)*yd(j,k);
dydt(j) += A(i,j)*yd(i,k);
k = k+1;
This eliminates the $O(n^3)$ matrix-vector ...
All-over this is a nicely structured code. The main problems are related to the Runge-Kutta solvers, where the first-order system was not uniformly applied to the computation.
What is obviously wrong can be found in the first lines of the solver file
def rk2_derivatives(edo, qk, pk, dt, bodies):
k1 = dt * edo(qk, pk, bodies)
k2 = dt * edo(qk + (dt * k1), ...
Based on the references you provided, I assume your initial conditions are random numbers (see second reference, Sect. IV)? If so, it's probably no surprise that you do not see higher-order convergence in time, at least not from the beginning. Your initial conditions are not smooth and it may take some time to get to a smooth solution/state. There, however, ...