> 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});

[Plot]

>