Satellite position computation using Runge-Kutta 4

Question

my issue is related to the Runge-Kutta 4 (RK4) method and the correct iteration steps required for the state vector of an orbiting satellite. The below code (in Python) describes the motion based on the description as per this link (http://www.navipedia.net/index.php/GLONASS_Satellite_Coordinates_Computation):

        if total_step_number != 0:   
            for i in range(1, total_step_number+1):                             
                #Calculate k1                
                k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
                k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
                k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
                 + ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
                 + XYZDDot[2]

                #Intermediate step to bridge k1 to k2
                XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
                XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
                XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
                XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
                XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
                XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
                radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))

                 ....

There is more to this code however I want to limit what I show to this for now. Basically, for those familiar with state vectors and using RK4, you can see that the position and velocity is updated in the intermediate step, but not the acceleration. My question is related to the calculation required in order to update too the acceleration. It would begin for example:

XYZDDot[0] = ...
XYZDDot[1] = ...
XYZDDot[2] = ...

...but what exactly comes after is not very clear. Any advice welcome!