Solve_ivp output, into orbital plot?

Question

I'm trying to plot an orbit of a moon around Jupiter using gravitational acceleration. I cannot seem to determine how to use the solve_ivp function appropriately. Something is just not clicking... I have created ODEs for a moon, related to Jupiter at the origin.

year = np.arange(0,31536000,3600)
G = 6.67408e-11
jupiter_xpos = 0
jupiter_ypos = 0
jupiter_vel = (0,0)
jupiter_mass = 1.89819e27
Io_orbit = 421700000
Io_xpos = -421700000
Io_ypos = 0
Io_xvel = 0
Io_yvel =  -1773400
Io_mass = 8.9319e22
Io = [Io_xpos,Io_xvel,Io_ypos,Io_yvel]


def ode(Moon,t_max):
    #Moon[0,1,2,3]=[x,v_x,y,v_y]
    output=[0,0,0,0]
    R = ((Moon[0]**2 + Moon[2]**2)**(3/2))
    output[0]=Moon[1]
    output[1]= -G*jupiter_mass*Moon[0]/R
    output[2]=Moon[3]
    output[3]= -G*jupiter_mass*Moon[2]/R
    return output

#This is where the problem is
sol= solve_ivp(ode,Io,year)


plt.plot(sol[:,0],sol[:,2],0,0,'ro')
plt.axis('equal')
plt.grid(True)
plt.show()

I'm hoping to achieve a 2D orbital plot like this...

and also track and plot each change in x and y position and velocity against time.

Answer 1

The documentation for solve_ivp shows that the parameters are

sol = solve_ivp(ode, [t0,tf], u0, tspan=year, atol = 100, rtol=1e-9)

where year=np.arange(t0,tf,hour). Then you find the solution values in sol.t for the repeated times and sol.y for the values.