ial conditions Cómo resolver ecuaciones de movimiento de primer orden de ODE con condiciones iniciales dadas.

He convertido una ecuación a una ecuación diferencial de primer orden y ahora quisiera resolver las ecuaciones de movimiento para múltiples períodos con condiciones dadas.

Las ecuaciones se resolverán con los siguientes valores iniciales:
x(0), y(0), vx(0), vy(0) = 3, 1, 2, 1.3 × π

Esta es la ecuación diferencial de primer orden, derivada de la ecuación de movimiento de la ley de gravitación universal de Newton para cuerpos celestes:

¿Cómo puedo resolver esto en Python? ¿Podría usar los métodos de Runge Kutta o de Kepler? Siento que algo no está funcionando bien.

import math
import matplotlib.pyplot as plt
import numpy as np

g = 6.674 * 10(-11)
M = 4(math.pi2)/g
x0 = 3 #posición x del centro o h
y0 = 1 #posición y del centro o k
vx0 = 2
vy0 = 1.3 math.pi

Intentando con el segundo orden

def RK2(y0, f, tlist):
    t = [tlist[0]]
    tf = tlist1
    dt = tlist[2]

y = [y0]

while t[-1] < tf:
    k1 = dt * f(y[-1],t[-1])
    k2 = dt * f(y[-1]+0.5*k1 , t[-1]+0.5*dt) ##Tomar la mitad del paso

    y.append( y[-1] + k2 )
    t.append(t[-1] + dt)
return(np.array(y), np.array(t) )