Whats wrong with my RK4 Orbital Propagator?

Question

I am working in a RK4 orbital propagator, and I have followed the basic structure of:

k1​=dt×f(tn​,yn​)

k2​=dt×f(tn​+dt/2​,yn​+k1/2​​)

k3​=dt×f(tn​+dt/2​,yn​+k2/2​​)

k4​=dt×f(tn​+dt,yn​+k3​)​

yn+1=yn+1/6(k1+2k2+2k3+k4)

tn+1=tn+dt

I have bassed this procces in a simpler one, created by me for RK1 implementation:

#NOTE: sat is an array with the structure [x,y,z,vx,vy,vz,ax,ay,az]
def pasosimple(sat,dt):
   
    #physics constants  
    G = 6.6742e-20
    m_s = 1.989e30
    mu=G*m_s

    #Position Change from previous values of velocity and acceleration
    dx=sat[3]*dt+sat[6]*(dt**2)/2
    dy=sat[4]*dt+sat[7]*(dt**2)/2
    dz=sat[5]*dt+sat[8]*(dt**2)/2

   #Velocity change due to previous acceleration
    dvx=sat[6]*dt
    dvy=sat[7]*dt
    dvz=sat[8]*dt

   #xyz update
    x=dx+sat[0]
    y=dy+sat[1]
    z=dz+sat[2]
   
    #With the xyz update, we calculate new accelerations
    ax=(-mu*(x)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    ay=(-mu*(y)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    az=(-mu*(z)/(np.sqrt((x)**2+(y)**2+(z)**2)**3))
    
    #Substraction to obtain the difference acceleration
    dax=ax-sat[6]
    day=ay-sat[7]
    daz=az-sat[8]
   

    dsat=np.array([dx,dy,dz,dvx,dvy,dvz,dax,day,daz])           
    sat=np.array([x,y,z,dvx+sat[3],dvy+sat[4],dvz+sat[5],ax,ay,az])
    return dsat,sat

This code works, as far as I know, and I have already tested it. Now, for imlementing the RK4, I'm doing this:

def rk4(sat,dt):
    d1,s1=pasosimple(sat,dt)            
    d2,s2=pasosimple(sat+d1/2,dt+0.5*dt)    
    d3,s3=pasosimple(sat+d2/2,dt+0.5*dt)
    d4,s4=pasosimple(sat+d3,dt+dt)     
    sat=sat+(1/6)*(d1+d2*2+d3*2+d4) 

    return sat

Which does not work. I was hoping that anyone could give some insight about what i am doing wrong. Thank you all.

EDIT:

I have tried this alternate configuration for pasosimple, following (or trying to follow) the comments below. The result did not work either. However, using only the new pasosimple code and not nesting it into rk4works nicely, which makes me think the problem is now in rk4.

def pasosimple(sat,dt):
    G = 6.6742e-20
    m_t=5.972e24
    m_s = 1.989e30
    m_m=6.39e23
    mu=G*m_s
    mu_t=G*m_t
    mu_m=G*m_m
    
    ax=(-mu*(sat[0])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3))
    ay=(-mu*(sat[1])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3)) 
    az=(-mu*(sat[2])/(np.sqrt((sat[0])**2+(sat[1])**2+(sat[2])**2)**3))
                        
    x_punto=np.array([sat[3],sat[4],sat[5],ax,ay,az]) 
    x=np.array([sat[0]+x_punto[0]*dt,sat[1]+x_punto[1]*dt,sat[2]+x_punto[2]*dt,sat[3]+x_punto[3]*dt,sat[4]+x_punto[4]*dt,sat[5]+x_punto[5]*dt])
    
    
    dsat=np.array([x_punto[0]*dt,x_punto[1]*dt,x_punto[2]*dt,x_punto[3]*dt,x_punto[4]*dt,x_punto[5]*dt,ax-sat[6],ay-sat[7],ay-sat[8]])
    sat=np.array([x[0],x[1],x[2],x[3],x[4],x[5],ax,ay,az])
    return dsat,sat