restart:f||1:= (x,y,z) -> 10*(y-x);
f||2:= (x,y,z) -> x*(28.0-z)-y;
f||3:= (x,y,z) -> x*y-8.0/3.0*z;with(plots);with(LinearAlgebra);N:=2000;
A := Matrix(N,3);
k1:=Matrix(1,3);
k2:=Matrix(1,3);
k3:=Matrix(1,3);
k4:=Matrix(1,3);
h:=0.01;
A[1,1]:=5;A[1,2]:=5;A[1,3]:=6;j:='j':
for j from 1 to N-1 do
x:=A[j,1]: y:=A[j,2]: z:=A[j,3]:
k:='k':
for k from 1 to 3 do
k1[1,k]:=h*f||k(x,y,z):
od:
k:='k':
for k from 1 to 3 do
k2[1,k]:=h*f||k(x+k1[1,1]/2,y+k1[1,2]/2,z+k1[1,3]/2):
od:
k:='k':
for k from 1 to 3 do
k3[1,k]:=h*f||k(x+k2[1,1]/2,y+k2[1,2]/2,z+k2[1,3]/2):
od:
k:='k':
for k from 1 to 3 do
k4[1,k]:=h*f||k(x+k3[1,1],y+k3[1,2],z+k3[1,3]):
od:
k:='k':
for k from 1 to 3 do
A[j+1,k]:= A[j,k]+(k1[1,k]+2*k2[1,k]+2*k3[1,k]+k4[1,k])/6.0:
od:
od:pointplot3d(A,connect=true,color = red );