clear all;

//Parameters
o = 10;
r = 28;
b = 8/3;

// x=x(1), y=x(2), z=x(3)
function dx = LORENZ(t,x)
  dx(1) = -o * x(1) + o * x(2);
  dx(2) = r * x(1) - x(2) - x(1) * x(3); 
  dx(3) = - b * x(3) + x(1) * x(2); 
endfunction

// The time range and step for calculation
Tin = 0;
Tout = 30;
Tstep = 0.01;
t = Tin:Tstep:Tout;

// The initial values
x10 = -4;
y10 = -6;
z10 = 15;
X10 = [x10; y10; z10];

X1 = ode(X10,0,t,LORENZ);

Tmax = Tout / Tstep
deltaT = 0.05
deltaI = deltaT / Tstep
Imax = Tmax - deltaI*2

K1 = zeros(1,Tmax);
K2 = zeros(1,Tmax);
K3 = zeros(1,Tmax);

for i = 1 : Imax
     K1(1, i) = [ X1(1, i) ] ;
     K2(1, i) = [ X1(1, i+deltaI) ] ;
     K3(1, i) = [ X1(1, i+deltaI*2) ] ;
end

scf(1);

param3d(K1(1,:), K2(1,:), K3(1,:))
xlabel('x1')
ylabel('x2')
zlabel('x3')

scf(2);
param3d(X1(1,:), X1(2,:), X1(3,:))