clear;

s = 10;
r = 28;
b = 8 / 3;

function dx = lorenz(t,x)
  dx(1) = - s * (x(1) - x(2));
  dx(2) = r * x(1) - x(2) - x(1) * x(3);
  dx(3) = x(1) * x(2) - b * x(3);
endfunction

T = linspace(0,30,3000);
x0 = 0;
y0 = 0.03;
z0 = 0;
X0 = [x0; y0; z0];

X = ode(X0,0,T,lorenz);
plot3d3(X(1,:),X(2,:),X(3,:));