I just need help getting this problem started in MATLAB, specifically what the ode_function should look like.

The satellite equation of motion can be described by: r = - mu r/r^3 + a_p where r = [x, y, z]^T is the position of the satellite, r = |r| middot mu = 3.98716708 times 10^5 km^3/s^2 is the gravitational parameter. a_p = [- 3 J_2 mu_2 R^2_ x/2r^5 (1 - 5z^2/r^2) -3J_2 mu R^2_ y/2r^(1 - 5z^2/r^2) - 3 J_2 mu R^2_ z/2r^5(3 - 5z^2/r^2)] is the perturbation due to the non spherical earth. R_ = 6378.137km is the radius of the earth. J_2 = 0.0010826267. Use different numerical integration methods to find the trajectory of the satellite assuming the initial conditions of r(0) = [-6035, -3480, 2500]^T km; r(0) = [-3, 437, 6.618, 2.523]^T km/s The integration time interval is set to [0 2.5]-hours. (1) Write a general Runge-Kutta integration function with this specified form: [tout, yout] = RK(ode_function, tspan, y0, h, rk) where the output and input parameters are defined as ode_function: user defined M-function in which the derivatives f are computed tspan: the vector of [t_0 t_r] giving the time interval for the solution, t_0 is the initial time and tr is the final time y0: column vector of initial values of the vector y tout: column vector of times at which y was evaluated yout: a matrix, each row of which contains r = [x, y, z] and r = [x, y, z] evaluated at the corresponding time in tout h: time step rk: rk can only take the values of 1, 2, or 4. rk = 1 indicates the first-order RK method or the Euler method, rk = 2 indicates the second-order RK method; rk = 4 indicates the fourth-order RK method. The satellite equation of motion can be described by: r = - mu r/r^3 + a_p where r = [x, y, z]^T is the position of the satellite, r = |r| middot mu = 3.98716708 times 10^5 km^3/s^2 is the gravitational parameter. a_p = [- 3 J_2 mu_2 R^2_ x/2r^5 (1 - 5z^2/r^2) -3J_2 mu R^2_ y/2r^(1 - 5z^2/r^2) - 3 J_2 mu R^2_ z/2r^5(3 - 5z^2/r^2)] is the perturbation due to the non spherical earth. R_ = 6378.137km is the radius of the earth. J_2 = 0.0010826267. Use different numerical integration methods to find the trajectory of the satellite assuming the initial conditions of r(0) = [-6035, -3480, 2500]^T km; r(0) = [-3, 437, 6.618, 2.523]^T km/s The integration time interval is set to [0 2.5]-hours. (1) Write a general Runge-Kutta integration function with this specified form: [tout, yout] = RK(ode_function, tspan, y0, h, rk) where the output and input parameters are defined as ode_function: user defined M-function in which the derivatives f are computed tspan: the vector of [t_0 t_r] giving the time interval for the solution, t_0 is the initial time and tr is the final time y0: column vector of initial values of the vector y tout: column vector of times at which y was evaluated yout: a matrix, each row of which contains r = [x, y, z] and r = [x, y, z] evaluated at the corresponding time in tout h: time step rk: rk can only take the values of 1, 2, or 4. rk = 1 indicates the first-order RK method or the Euler method, rk = 2 indicates the second-order RK method; rk = 4 indicates the fourth-order RK method