x,y,z,w,t=var('x y z w t')
a=10
b=28
c=8/3
init_x=2
init_y=3
init_z=4
init_w=5
p=desolve_system_rk4([a*(y-x)+sin(w*t),x*(b-z)-y+sin(w*t),x*y-c*z+sin(w*t),sin(3*w*t)],[x,y,z,w],ics=[0,init_x,init_y,init_z,init_w],ivar=t,end_points=20)
q1=[[i,j] for i,j,k,1,2 in p]
lp1=list_plot(q1,plotjoined=true)
q2=[[i,k] for i,j,k,1,2 in p]
lp2=list_plot(q2,plotjoined=true, color='red')
q3=[[i,1] for i,j,k,1,2 in p]
lp3=list_plot(q3,plotjoined=true, color='green')
q4=[[i,2] for i,j,k,1,2 in p]
lp4=list_plot(q4,plotjoined=true, color='blue')
lp1+lp2+lp3
