function dy = pam(t,y,p)
% p = [s0, sigma, mu0, Smax, beta, theta, I, e0, c];
% y = [S, N, P];
% Example:
%   global P;
%   P = [1 1 0.0001 200 0.0004 1 0.01 5 0.04];
%   y0 = [0 0 0];
%   t = 0:10:1000;
%   opts = odeset('RelTol',1e-3,'AbsTol',1e-6);
%   [t,y] = ode45('pam', t, y0, opts);
%   plot(t,y(:,1),'-');
%   plot(y(:,3),y(:,1),'-');
global P;
dy = zeros(3,1);    % a column vector
dy(1) = P(1)*P(2) + (P(4)*P(5) - P(1)*P(2))*y(2) - (P(5) + P(3) + P(6)*y(3)*exp(-P(9)*y(1)))*y(1);
dy(2) = P(2) - (P(2) + P(3) + P(6)*y(3)*exp(-P(9)*y(1)))*y(2);
if y(1) == 0
   dy(3) = 0;
else
   dy(3) = P(7) - (P(8)/P(6))*(1/y(1))*exp(P(9)*y(1))*y(3);
end