   function [Enter1,Exit1,Enter2,Exit2,tOut,yOut] = LorenzZones(y0,tFinal)
% Integrates the IVP
%
%      [ y1'(t) ]      [  beta     0        y2(t) ]  [ y1(t) ]
%      [ y2'(t) ]  =   [   0     -sigma    sigma  ]  [ y2(t) ]    y(0) = y0
%      [ y3'(t) ]      [ -y2(t)  rho         -1   ]  [ y3(t) ]
%
% from t = 0 to t= tFinal. The parameter values are...

  sigma = 10;
  rho = 28;
  beta = 8/3;
  
% Think of y(t) as the location of a comet at time t.
% Uses ode45 with Relerr = .000001. 
%
% Enter1 a column vector (possibly empty) indicating when the comet enters Zone 1
% Exit1  a column vector (possibly empty) indicating when the comet leaves Zone 1
% Enter2 a column vector (possibly empty) indicating when the comet enters Zone 2
% Exit2  a column vector (possibly empty) indicating when the comet leaves Zone 2
%
% the 3-column matrix yOut houses the ode45 solution snapshops associated with 
% the times in the vector tOut.


% The two attractors: 
eta = sqrt(beta*(rho-1));
u = [rho-1;eta;eta];           % Attractor 1
v = [rho-1;-eta;-eta];         % Attractor 2

Enter1 = [];
Exit1  = [];
Enter2 = [];
Exit2  = [];
opts = odeset('RelTol',.000001,'events',@(t,y) LorenzEvents(t,y,u,v));
[tOut,yOut,tE,yE,iE]= ode45(@(t,y) Lorenz(t,y,beta,sigma,rho),[0 tFinal],y0,opts);
Enter1  = tE(find(iE==1));
Exit1   = tE(find(iE==2));
Enter2  = tE(find(iE==3));
Exit2   = tE(find(iE==4));

  
function yp = Lorenz(t,y,beta,sigma,rho)
yp = [-beta*y(1)+y(2)*y(3)    ;...
      -sigma*y(2)+sigma*y(3)  ; ...
      -y(2)*y(1)+rho*y(2)-y(3)...
      ];
  
  
 function [gstop,isterminal,direction] = LorenzEvents(t,y,u,v)
% Event function for turning off ode45 after (approximately) one
% revolution. We do this when dx/dt is zero and decreasing.

% d is positive  if the comet is closer to u = Attractor 1...
d =  sqrt((v(1)-y(1))^2 + (v(2)-y(2))^2 + (v(3)-y(3))^2) - ...
     sqrt((u(1)-y(1))^2 + (u(2)-y(2))^2 + (u(3)-y(3))^2);
% c is positive if angle yuv is acute
c = (y-u)'*(v-u);

% Comet is in Zone 1 if d>0
% Comet is in Zone 2 if c<0
% Event definitions:
%    1  Comet entering Zone 1
%    2  Comet leaving  Zone 1
%    3  Comet entering Zone 2
%    4  Comet leaving  Zone 2


gstop      = [  d d  c c];
direction  = [ 1 -1 -1 1];
isterminal = [ 0 0 0 0];
