% function to specify the set of ODEs to be integrated by ode45
%
function dy = gause(t,y)               %gause will be the name of the 
                                       %function to be called by 1st argument
                                       %of ode45

dy = zeros(3,1);                       %this says that dy is a vector of 3
                                       %components and initialize them to 
				       %zero

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin specifying set of ODEs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Suggested parameters to keep fixed
alpha=1.5;  % Crowding advantage in carrying capacity of N_1 over N_2
b=0.01;     % Predation rate of N_1
c=0.5;      % Typical per-predator increase when eating one prey at a rate b (or b-eps)
d=1;	    % Mortality rate of predators
r=1;        % Reproductive rate of prey

% Modify these
eps=0.009;  % Predator avoidance advantage of N_2 versus N_1
k=925;      % Carrying capacity of N_2

dy(1)=y(1)*(r-r/k*y(1)-r/k*y(2)-b*y(3));
dy(2)=y(2)*(r-r/k*alpha*y(1)-r/k*y(2)-(b-eps)*y(3));
dy(3)=y(3)*(c*b*y(1)+c*(b-eps)*y(2)-d);
