function [T,Y]=lorenz_ode45(xo,yo,zo,tf)
%
%Integrates the Lorenz systems
%
%  Usage: [T,Y]=lorenz_ode45(x0,y0,z0,tf); 
%
%  Example:
%             [T,Y]=lorenz_ode45(8,3,33,50);


%options= odeset('Reltol',1.e-5);
options= odeset('Reltol',1.e-3,'OutputFcn','odephas3','Refine',4);

 figure(1); clf;
 IC=[xo;yo;zo];
 tspan=[0 tf];
 
 [T,Y]=ode45(@lorenz_eqns,tspan,IC,options);
 
 % ploting the orbit in 3D
 figure(2); clf
 plot3(Y(:,1),Y(:,2),Y(:,3))
 hold on
 xlabel('x')
 ylabel('y')
 zlabel('z')
 title('Lorenz Attractor')
 rotate3d
 view([60 45]);
 axis tight
 hold off

 % ploting the time series for each coordinate
 figure(3); clf
 hold on
 xlabel('t')
 ylabel('x!y!z')
 plot(T,Y(:,1),'r-',T,Y(:,2),'b-',T,Y(:,3),'g-')
 hold off
 axis tight

 return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Differential equations defining the Lorenz system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=lorenz_eqns(t,vec)

 sigma=10.0;
 rho=28.0;
 beta=8.0/3.0;

 x=vec(1);
 y=vec(2);
 z=vec(3);

 F(1)= sigma* (y -x);
 F(2)= rho* x - y - x*z;
 F(3)= -beta* z + x*y;

 F=transpose(F);
 return;

