Orbit spirals using 4th order Yoshida integration

Question

I am attempting to use the 4th order Yoshida integration technique to model the orbit of satellites in circular orbits around the Earth.

However, the orbits I achieve spiral away quite quickly. The code for a Moon like satellite is below. Interestingly, the particles behaved when I use Euler method, however, I wanted to try a more accurate method. The issue could then be within how I have implemented the algorithm itself.

I have tried using the gravitational parameter rather then computing G*M, but this did not help. I also reduced the time-step, messed around with units, printed and checked values for various integration steps etc., but could not find anything. Is this the correct use of this algorithm?

G = 6.674e-20 # km^3 kg^-1 s^-2
day = 60.0 * 60.0 * 24.0 # length of day
dt = day / 10.0
M = 5.972e24 # kg

N = 1
delta = np.random.random(1) * 2.0 * np.pi / N
angles = np.linspace(0.0, 2.0 * np.pi, N) + delta
rad = np.random.uniform(low = 384e3, high = 384e3, size = (N))

x, y = rad * np.cos(angles), ringrad * np.sin(angles)
vx, vy = np.sqrt(G*M / rad) * -np.sin(angles), np.sqrt(G*M / rad) * np.cos(angles)

def update(frame):
    global x, y, vx, vy, dt, day

    positions.set_data(x, y)

    # coefficients 
    q = 2**(1/3)
    w1 = 1 / (2 - q)
    w0 = -q * w1
    d1 = w1
    d3 = w1
    d2 = w0
    c1 = w1 / 2
    c2 = (w0 + w1) / 2
    c3 = c2
    c4 = c1

    # Step 1

    x1 = x + c1*vx*dt
    y1 = y + c1*vy*dt
    dist1 = np.hypot(x1, y1)
    acc1 = -(G*M) / (dist1**2.0)
    dx1 = x1 - x
    dy1 = y1 - y
    accx1 = (acc1*dx1)/(x1)
    accy1 = (acc1*dy1)/(y1)
    vx1 = vx + d1*accx1*dt
    vy1 = vy + d1*accy1*dt

    # Step 2

    x2 = x1 + c2*vx1*dt
    y2 = y1 + c2*vy1*dt
    dist2 = np.hypot(x2, y2)
    acc2 = -(G*M) / (dist2**2.0)
    dx2 = x2 - x1
    dy2 = y2 - y1
    accx2 = (acc2*dx2)/(x2)
    accy2 = (acc2*dy2)/(y2)
    vx2 = vx1 + d2*accx2*dt
    vy2 = vy1 + d2*accy2*dt

    # Step 3

    x3 = x2 + c3*vx2*dt
    y3 = y2 + c3*vy2*dt
    dist3 = np.hypot(x3, y3)
    acc3 = -(G*M) / (dist3**2.0)
    dx3 = x3 - x2
    dy3 = y3 - y2
    accx3 = (acc3*dx3)/(x3)
    accy3 = (acc3*dy3)/(y3)
    vx3 = vx2 + d3*accx3*dt
    vy3 = vy2 + d3*accy3*dt

    # Full step

    x = x3 + c4*vx3*dt
    y = y3 + c4*vy3*dt
    vx = vx3
    vy = vy3

    return positions