% predpre.m - discrete time predator prey example
% RMM, 6 Oct 03

aminit;					% initialize MATALB

%
% Predator prey system
%

% Set up the initial state
clear H L year
% H(1) = 20; L(1) = 35;
H(1) = 10; L(1) = 10;

% For simplicity, keep track of the year as well
year(1) = 1845;

% Set up parameters
br = 0.6; df = 0.7; a = 0.014;
nperiods = 365;
duration=90;

% Iterate the model
for k = 1:duration*nperiods
  b = br;					% constant food supply
% b = br*(1+0.5*sin(2*pi*k/(4*nperiods)));	% varying food supply (try it!)
  H(k+1) = H(k) + (b*H(k) - a*L(k)*H(k))/nperiods;
  L(k+1) = L(k) + (a*L(k)*H(k) - df*L(k))/nperiods;
  year(k+1) = year(k) + 1/nperiods;

  if (mod(k, nperiods) == 1)
    % Store the annual population
    Ha((k-1)/nperiods + 1) = H(k);
    La((k-1)/nperiods + 1) = L(k);
  end;
end;

% Store the final population
Ha(duration) = H(duration*nperiods+1);
La(duration) = L(duration*nperiods+1);

% Plot the populations of rabbits and foxes versus time
subplot(211)
h = plot(1845 + [1:duration], Ha, '.-', 1845 + [1:duration], La, '.--');
set(h, 'Linewidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize*3);
amaxis([1845 1925 0 250]);

xlabel('Year');
ylabel('Population');

% Now reset the parameters to look like we want
daspect([12 125 1]);
lgh = legend(gca, 'hares', 'lynxes', 'Location', 'NorthEast', ...
  'Orientation', 'Horizontal');
legend(lgh, 'boxoff');

amprint('predprey-discrete.eps');

