from sage.calculus.desolvers import desolve_system_rk4 
x,y,vx,vy,t=var('x y vx vy t')
def r(x,y):
    return  sqrt(x**2+y**2)
def fx(x,y):
    return -x/(r(x,y)**3)
    
def fy(x,y):
    return  -y/(r(x,y)**3)
    
       
P=desolve_system_rk4([vx,fx(x,y),vy, fy(x,y)],[x,vx,y,vy],ics=[0,1,0,1,1],step=0.01,ivar=t,end_points=30)
P1=desolve_system_rk4([vx,fx(x,y),vy, fy(x,y)],[x,vx,y,vy],ics=[0,1,0,1,0.8],step=0.01,ivar=t,end_points=30)
P2=desolve_system_rk4([vx,fx(x,y),vy, fy(x,y)],[x,vx,y,vy],ics=[0,1,0,1,-1.2],step=0.01,ivar=t,end_points=90)
P3=desolve_system_rk4([vx,fx(x,y),vy, fy(x,y)],[x,vx,y,vy],ics=[0,1,0,1,1.2],step=0.01,ivar=t,end_points=90)
Q=[[x,y] for t,x,vx,y,vy in P]
Q1=[[x,y] for t,x,vx,y,vy in P1]
Q2=[[x,y] for t,x,vx,y,vy in P2]
Q3=[[x,y] for t,x,vx,y,vy in P3]
R=list_plot(Q,plotjoined=True,rgbcolor=(1,0,0))
R1=list_plot(Q1,plotjoined=True,rgbcolor=(0,1,0))
R2=list_plot(Q2,plotjoined=True,rgbcolor=(0,0,1))
R3=list_plot(Q3,plotjoined=True,rgbcolor=(0,0,1))
R+R1+R2+R3