solving delay differential equations, Heun's method with circular buffer

Question

Considering the following toy example for solving DDE using Heun's method, I need to verify if Heun's method is functioning correctly.

I have visually compared Heun's method and Euler's method, as well as the value of y at corresponding times. Additionally, it's important to note that the time step in Heun's method is one order of magnitude larger. The error, in comparison to Euler's method, need to be of the order O(dt)^2.

In the function fn, y represents the buffer, while yn denotes the latest coordinates.

import numpy as np 
# from numba import njit
import matplotlib.pyplot as plt
from numpy.random import rand, randn

# @njit
def fn(yn, y, D):
    f0 = -yn[0] * y[1, -D[0]] + (y[1, -D[1]])
    f1 = yn[0] * y[1, -D[0]] - yn[1]
    f2 = yn[1] - y[1, -D[1]]

    return np.array([f0, f1, f2])

#@njit
def euler_step(y, D, dt):    

    yn = y[:, -1]
    k = fn(yn, y, D)
    new_y = yn + dt * k
    return new_y

# @njit
def heun_step(y, D, dt):
    yn = y[:, -1]
    k1 = fn(yn, y, D)
    y1 = yn + dt * k1
    k2 = fn(y1, y, D)
    new_y = yn + 0.5 * dt * (k1 + k2)
    return new_y
    
# @jit(nopython=False)
def dde_solve(L, y0, t, dt, seed=2, decimate=10, method='euler'):
    
    np.random.seed(seed)
    nstep = len(t)
    D = np.round(L/dt).astype(np.int64)
    md = np.max(D)
    y = np.tile(y0, (md, 1)).T
    
    store = np.zeros((len(y0), nstep//decimate))
    for it in range(nstep):
        if method == 'euler':
            new_y = euler_step(y, D, dt)
        elif method == 'heun':
            new_y = heun_step(y, D, dt)
        
        y[:, :-1] = y[:, 1:]
        y[:, -1] = new_y
        # store
        if it % decimate == 0:
            store[:, it//decimate] = new_y
        
    return store
    
seed = 0
tend = 20
decimate= 10
dt = 0.01
L = np.array([1.0, 10.0])
t = np.arange(0.0, tend, dt)
y0 = np.array([5.0, 0.1, 1.0])
    
"euler:"
y = dde_solve(L, y0, t, dt, seed, decimate, 'euler')
plt.plot(y[0,:], y[1,:], 'r', label='euler', alpha=0.5)
i = np.where(np.abs(t[::decimate] - 10) < 1e-6)[0][0]
print("euler ", y[:, i], t[::decimate][i])

"heun:"
dt = 0.1
t = np.arange(0.0, tend, dt)
decimate = 1
y = dde_solve(L, y0, t, dt, seed, decimate, 'heun')
plt.plot(y[0,:], y[1,:], 'b', label='heun', alpha=0.5, ls="-.")
i = np.where(np.abs(t[::decimate] - 10) < 1e-6)[0][0]
print("heun ", y[:, i], t[::decimate][i])

plt.legend()
plt.show()


# Comparing the results for 3 different dt, in heun dt is 10 times larger.
# euler  [0.33621002 0.03793458 5.7258554 ] 10.0  0.01
# heun   [0.35603672 0.03547178 5.7084915 ] 10.0  0.1

# euler  [0.3331855  0.03905495 5.72775955] 10.0  0.001
# heun   [0.33484529 0.03875237 5.72640234] 10.0  0.01

# euler  [0.33288846 0.03916828 5.72794326] 10.0  0.0001
# heun   [0.33305097 0.03913745 5.72781157] 10.0  0.001

Observing the comparisons, I'm unsure whether the error order and time steps are compatible.