% Predator prey example
%
% The code implementing this is in the file predprey.m.  The parameter
% values for the plots generate in this lecture are
%
%

aminit;			% initialize MATLAB

alpha = 12.5;		% multiplier for hare population
beta = 25;		% multipler for lynx population
gamma = 0.8;		% multiplier for time scale

r = 2 * gamma;			% birth rate of hares
d = 0.7 * gamma;		% death rate for lynxes
b = 0.3 * beta / alpha;		% birth coefficient for lynxes
k = 10 * alpha;			% carrying capacity for hares
a = 8 * gamma * alpha / beta;	% interaction coefficient
c = 4 * alpha;			% saturation coefficient

x0 = [25, 20];
Tmax = 90;

% Generate a simulation of the system
figure(1); clf; subplot(211);
[t, x] = ode45('predprey', [0 Tmax], x0, [], r, d, b, k, a, c);
xy = plot(t, x(:,1), '-', t, x(:,2), '--');
set(xy, 'LineWidth', AM_data_linewidth);
amaxis;
xlabel('Time {\itt} [years]');
ylabel( 'Population', 'Rotation', 90);

% Pick out the periodic solution
Hf = x(length(t),1);
Lf = x(length(t),2);

lgh = legend('Hare', 'Lynx', 'Location', 'Northwest'); legend(lgh, 'boxoff');
amaxis([0 70 0 100])

daspect(daspect * diag([1 1.5 1]));

amprint('predprey-sim.eps');

% Generate the phase portrait
figure(2); clf; subplot(221); hold on;

% Generate a vector field
amphaseplot('predprey', [0 60 7], [0 50 6], 0.2, [], ...
 0, r, d, b, k, a, c); 
amphaseplot('predprey', [0 60 7], [60 100 4], 0.06, [], ...
 0, r, d, b, k, a, c);
amphaseplot('predprey', [70 120 6], [0 50 6], 0.075, [], ...
 0, r, d, b, k, a, c);
amphaseplot('predprey', [70 120 6], [60 100 4], 0.05, [], ...
 0, r, d, b, k, a, c);

% Plot the limit cycle
[ah, lh] = amphaseplot('predprey', 'auto', 19, 0, ...
  [Hf,Lf], 11, r, d, b, k, a, c);
set(lh, 'LineWidth', AM_thick_linewidth, 'Color', [1 0 0]);

% Plot some outside trajectories
[ah, lh] = amphaseplot('predprey', 'auto', 19, 0, ...
  [120,32; 120,60], ...
  18, r, d, b, k, a, c);
set(lh, 'LineWidth', AM_data_linewidth);

% Plot some inside trajectories
[ah, lh] = amphaseplot('predprey', 'auto', 19, 0, ...
  [19, 30], ...
  75, r, d, b, k, a, c);
set(lh, 'LineWidth', AM_ref_linewidth);

amaxis([-1 120 -1 100]);
xlabel('Hares'); ylabel('Lynxes', 'Rotation', 90);

% Compute the location of the equilibrium point and plot it
Xeq = fsolve('predprey_ss', [20, 30], [], r, d, b, k, a, c);

hold on; plot([0, Xeq(1)], [0, Xeq(2)], 'k.', ...
  'MarkerSize', AM_data_markersize*3);
amprint('predprey-pp.eps');
