from sage.calculus.desolvers import desolve_system_rk4 
x,y,t=var('x y t')
P=desolve_system_rk4([y,-sin(x)],[x,y],ics=[0,0.1,2.001],step=0.01,ivar=t,end_points=20)
P2=desolve_system_rk4([y,-sin(x)],[x,y],ics=[0,0.1,1.99],step=0.01,ivar=t,end_points=20)
P1=desolve_system_rk4([y,-x],[x,y],ics=[0,0.1,2.001],step=0.01,ivar=t,end_points=20)
Q=[ [j,k] for i,j,k in P]
Q1=[ [j,k] for i,j,k in P1]
Q2=[ [j,k] for i,j,k in P2]
R1=list_plot(Q1,plotjoined=True,rgbcolor=(1.0,0,0),aspect_ratio=1)
R=list_plot(Q,plotjoined=True,rgbcolor=(0.0,1.0,0),aspect_ratio=1)
R2=list_plot(Q2,plotjoined=True,rgbcolor=(0.0,0.1,0),aspect_ratio=1)
R1+R+R2