function ShowMultipleEvent()
% Use the rabbit-fox simulation to illustrate the use of
% multiple event functions

close all

% The population densities for rabbits and foxes at time t=0...
r0 = 300;
f0 = 150;
y0 = [r0;f0];
t0 = 0;
tFinal = 20;
alpha = .01;

% Integrate across the [t0,tFinal] and highlight all the
% local mins for the rabbits
options = odeset('events',@(t,y) g(t,y,alpha));
[tVals,yVals,tE,yE,iE] = ode23(@(t,y) MyF(t,y,alpha),[t0 tFinal],y0,options);
figure
plot(tVals,yVals(:,1),'k','linewidth',2)
hold on
plot(tVals,yVals(:,2),'r','linewidth',2)
%axis([0 20 0 500])
legend('Rabbits','Foxes')

% Highlight the rabbit maximums (Event 1)...
idx1 = find(iE==1);
plot(tE(idx1),yE(idx1,1),'.c','markersize',25)

% Highlight the rabbit minimums (Event 2)...
idx2 = find(iE==2);
plot(tE(idx2),yE(idx2,1),'.m','markersize',25)

% Highlight the fox maximums (Event 3)...
idx3 = find(iE==3);
plot(tE(idx3),yE(idx3,2),'.b','markersize',25)

% Highlight the fox minimums (Event 4)...
idx4 = find(iE==4);
plot(tE(idx4),yE(idx4,2),'.g','markersize',25)
hold off
shg

% Rabbits go extinct..
if iE(length(iE))== 5
xlabel('Rabbits Go Extinct','fontsize',24)
end

% Foxes go extinct..
if iE(length(iE))== 6
xlabel('Foxes Go Extinct','fontsize',24)
end

function yp = MyF(t,y,alpha)
% y1 = rabbits, y2 = foxes
yp = [2*y(1)-alpha*y(1)*y(2) ; -y(2)+alpha*y(1)*y(2)];

function [gStop,isterminal,direction] = g(t,y,alpha)
% Event function for rabbit-fox simulation.
% Six events; how to detect them, what to do if they happen:
%
% Event What Trigger Keep Integrating?
% 1 Local max for the rabbits y1'(t)=0, y1'(t) increasing yes
% 2 Local min for the rabbits y1'(t)=0, y1'(t) decreasing yes
% 3 Local max for the foxes y2'(t)=0, y2'(t) increasing yes
% 4 Local min for the foxes y2'(t)=0, y2'(t) decreasing yes
% 5 Rabbits almost extinct y1(t)=1 decreasing no
% 6 Foxes almost extinct y2(t)=0 decreasing no

yp = MyF(t,y,alpha);
% The vector of triggers...
gStop = [yp(1); yp(1); yp(2); yp(2); y(1)-1; y(2)-1];
% The vector of tell-tale directions...
direction = [-1; 1; -1; 1;-1; -1];
% The "do we halt the integration" vector...
isterminal = [0; 0; 0 ; 0; 1;1];