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=60)
Q=[[x,y] for t,x,vx,y,vy in P]
R=list_plot(Q,plotjoined=True,rgbcolor=(0,0,1))
movingpoints=[R+point( (Q[47*k][0],Q[47*k][1]), color='red',size=80 ) for k in range(0,79)]
a=animate(movingpoints,xmin=-5,xmax=5,ymin=-5,ymax=5)
a.show()