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

 params=[10,28,8/3]; % params=[sigma,rho,beta]
 tspan=[0 tf];

 options= odeset('Reltol',1.e-5);
 %options= odeset('Reltol',1.e-3,'OutputFcn','odephas3','Refine',4); figure(1); clf;
 
 [T,Y]=ode45(@(t,y) lorenz_RHS(t,y,params),tspan,IC,options);
 
 % ploting the orbit in 3D
 figure(2); clf
 plot3(Y(:,1),Y(:,2),Y(:,3))
 xlabel('x')
 ylabel('y')
 zlabel('z')
 title('Lorenz Attractor')
 view([60 45]);
 axis tight

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

 return;


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

 sigma=params(1);
 rho  =params(2);
 beta =params(3);

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

 F1= sigma* (y -x);
 F2= rho* x - y - x*z;
 F3= -beta* z + x*y;

 F=[F1;F2;F3];
 return;

