>

restart;

>

>	f1:=(t,x,y,z)->10*y-10*x;

>	f2:=(t,x,y,z)->28*x-y-x*z;

>	f3:=(t,x,y,z)->-8/3*z+x*y;

>	t0:=0:x0:=5:y0:=5:z0:=5:

>	XY:=[x0,y0]:XYZ:=[x0,y0,z0]:

>	h:=0.01:n:=2000:t:=t0:x:=x0:y:=y0:z:=z0:

>	for i from 1 to n do

>	a:=f1(t,x,y,z):b:=f2(t,x,y,z):c:=f3(t,x,y,z):

>	u:=x+h*a:v:=y+h*b:w:=z+h*c:

>	x:=x+h*(a+f1(t+h,u,v,w))/2:y:=y+h*(b+f2(t+h,u,v,w))/2:z:=z+h*(c+f3(t+h,u,v,w))/2:

>	XY:=(XY,[x,y]):XYZ:=(XYZ,[x,y,z]):

>

t:=t+h:

>

od:

>	with(plots):

>	pointplot3d({XYZ});

>