function ShowEvent()
% Use the rabbit-fox simulation to illustrate the use of 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;

% Integrate across the [t0,tFinal] and display...

[tVals,yVals] = ode23(@MyF,[t0 tFinal],y0);
plot(tVals,yVals(:,1),'k','linewidth',2)
hold on
plot(tVals,yVals(:,2),'r','linewidth',2)
legend('Rabbits','Foxes')
hold off

% Integrate across the [t0,tFinal] and stop at the first local
% max in the rabbit density that occurs after t=8

options = odeset('events',@g);
[tVals,yVals] = ode23(@MyF,[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')
hold off

% Integrate across the [t0,tFinal] and highlight all the
% local mins for the rabbits
options = odeset('events',@g1);
[tVals,yVals,tE,yE] = ode23(@MyF,[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')
plot(tE,yE(:,1),'.c','markersize',25)
hold off

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

function [gStop,isterminal,direction] = g(t,y)
% Event function for turning off ode23 when the
% rabbit population is at a local max and t>=8
yp = MyF(t,y);
if t>8
gStop = yp(1); % This is the rate of change function for rabbits
else
gStop = 1; % Ensures that gStop is not zero before t=8.
end
isterminal = 1; % Stop ode23 when the event happens
direction = -1; % At a local max, the yp(1) function should be decreasing

function [gStop,isterminal,direction] = g1(t,y)
% Event function that looks for local mins for the rabbits
yp = MyF(t,y);
gStop = yp(1);
isterminal = 0; % Make a note of the event, but keepintegrating
direction = 1; % At alocal min, yp(1) should be increasing