%Lorenz equations
% Inputs: interval [a,b], initial vector y0
%
% Example usage: lorenz1([0 10],[1 1 10]);
function lorenz1(int,y0)
tol=odeset('Reltol', 1e-6);
[t y]=ode45(@ydot,int,y0,tol);
h=1e-10;
dt=0.01;
[t1 y1]=ode45(@ydot,int,y0+[h 0 0],tol);
%comet3(y(:,1),y(:,2),y(:,3));
%comet(y(:,1), y(:,2))
figure(1)
plot(y(:,1), y(:,2))
xlabel('y(1)')
ylabel('y(2)')
t0=int(1):dt:int(2);
yy=interp1(t,y,t0);
yy1=interp1(t1,y1,t0);
figure(2)
plot(t0,yy(:,1)-yy1(:,1))
xlabel('t'); ylabel('y1(1)-y(1)');
figure(3)
hold on
axis([-30 30 -30 30]);
for i=1:size(t0,2)
    plot(yy(i,1),yy(i,2),'-')
    pause(1e-100);
end
% figure(3)
% plot(y(:,2))
% figure(4)
% plot(y(:,3))

%plot(y)
%plot::Curve3d(y, t=int(1)..int(2))

function z=ydot(t,y)
% Lorenz equations
s=10; r=28; b=8/3;
z=zeros(3,1);
z(1)=-s*y(1)+s*y(2);
z(2)=-y(1)*y(3)+r*y(1)-y(2);
z(3)=y(1)*y(2)-b*y(3);