Resuelve las ecuaciones de movimiento para la primera ODE utilizando scipy.

Me gustaría resolver ecuaciones ODE de primer orden de movimiento utilizando la función “solve_ivp” de la librería “scipy”. Puedo ver que estoy haciendo algo mal ya que esto debería ser una elipse, pero solo estoy graficando cuatro puntos. ¿Puedes identificar el error?

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

gim = 4*(math.pi**2)
x0 = 1 #posición en x del centro o h
y0 = 0 #posición en y del centro o k
vx0 = 0 #posición vx
vy0 = 1.1*2*math.pi #posición vy
initial = [x0, y0, vx0, vy0] #estado inicial del sistema
time = np.arange(0, 1000, 0.01) #período

def motion(t, Z): 
  dx = Z[2] #vx
  dy = Z[3] #vy
  dvx = -gim/((Z[0]**2)+(Z[1]**2))**(3/2) * Z[0]
  dvy = -gim/((Z[0]**2)+(Z[1]**2))**(3/2) * Z[1]
  return [dx, dy, dvx, dvy]

sol = scipy.integrate.solve_ivp(motion, t_span=time, y0=initial, method='RK45')
plt.plot(sol.y[0],sol.y[1],"x", label="Solución Scipy RK45")
plt.show()

La trama que obtengo:

La trama que busco:

“`