function [P,Z,t] = pz(Po,Zo,T,method,dt)
% PZ(Po,Zo,T,method,dt)
%
%	Po	- initial value of P (at time=0)
%	Zo	- initial value of Z (at time=0)
%	T	- duration (time)
%	method	- integration method ('heun','rk4','forward',...)
%	dt	- time step
%
% Integrates a P-Z model with parameterized N using different time-stepping
% methods:
%
%                 N = Nt - P - Z
%             d_t P = u*P*N/(N+No) - g*Z*P/f(P+Po)
%             d_t Z = a*g*Z*P/f(P+Po) - d*Z
%
% Model parameters are set in the sub-function (in this script) called model()
% with values u=0.03, No=0.1, Nt=5, g=0.2, a=0.4, d=0.08. f(P)=1 is used.
% Model and parameters courtesy of Glenn Flierl and Stephanie Dutkiewicz.
%
% Possible values for 'method' are:
%	'forward'	- Forward method
%	'matsuno'	- Matsuno method (forward/backward)
%	'heun'		- Huen (forward/trapezoidal)
%	'rk2'		- Runge-Kutta 2 (mid-point method)
%	'rk4'		- Runge-Kutta 4
%	'ab2'		- Adams-Bashforth 2 (extrapolates dPdt,dZdt in time)
%	'ab2s'		- Adams-Bashforth 2 (extrapolates P,Z in time)
%	'ab3'		- Adams-Bashforth 3 (extrapolates dPdt,dZdt in time)
%	'ab3s'		- Adams-Bashforth 3 (extrapolates P,Z in time)
%	'qab2'		- (Quasi-) Adams-Bashforth 2 (extrap. dPdt,dZdt in time)
%	'qab2s'		- (Quasi-) Adams-Bashforth 2 (extrapolates P,Z in time)
%
% e.g.
%
% >> [P,Z,t]=pz(2,.5,500,'rk4',5);
% >> clf;plot(t,[P Z]);legend('P','Z');
% >> axes('position',[.21 .65 .25 .25]);plot(P,Z);xlabel('P');ylabel('Z');

% Number of time steps
nt=max(1,round(T/dt));

P=zeros(nt,1);	% Declares memory (faster to do this here)
Z=zeros(nt,1);	% Declares memory (faster to do this here)

% Initial conditions
P(1)=Po;
Z(1)=Zo;

for n=1:nt-0;
 nl=max(1,n-2:n);
 [P(n+1),Z(n+1)] = integrator( P(nl), Z(nl), dt, method );
end
t=dt*(0:nt-0)';

% -----------------------------------------------------------------------------

function [dPdt,dZdt] = model(P,Z,method)
% returns rates of change, d/dt P and d/dt Z as function of P and Z
%
% P and Z each contain two elements corresponding to time levels (n-1) and (n).

% Model parameters
Nt=5;		% Total nutrient
No=0.1;		% Threshold nutrient
u=0.03;		% Uptake factor (or growth rate of P)
g=0.2;		% Grazing factor
a=0.4;		% Grazing conversion factor (efficiency)
d=0.08;		% Death rate (of Z)
ff=1;

%u=0.01;	% Uptake factor (or growth rate of P)
%g=0.1;		% Grazing factor
%a=1.0;		% Grazing conversion factor (efficiency)
%Po=0.5;	% Threshold phytoplankton
%ff=P+Po;

N=Nt-P-Z;
dPdt = u*P*N/(N+No) - g*Z*P/ff;
dZdt = a*g*Z*P/ff - d*Z;

% -----------------------------------------------------------------------------

function [P,Z] = integrator(P,Z,dt,method)

Pn_2=P(end-2);
Pn_1=P(end-1);
P=P(end);
Zn_2=Z(end-2);
Zn_1=Z(end-1);
Z=Z(end);

switch method
 case 'forward',
  [dPdt,dZdt]=model(P,Z);
  P=P+dt*dPdt; Z=Z+dt*dZdt;
 case 'matsuno',
  [dPdt1,dZdt1]=model(P,Z); P1=P+dt*dPdt1; Z1=Z+dt*dZdt1;
  [dPdt2,dZdt2]=model(P1,Z1);
  P=P+dt*dPdt2; Z=Z+dt*dZdt2;
 case 'heun',
  [dPdt1,dZdt1]=model(P,Z); P1=P+dt*dPdt1; Z1=Z+dt*dZdt1;
  [dPdt2,dZdt2]=model(P1,Z1);
  P=P+dt*(dPdt1+dPdt2)/2; Z=Z+dt*(dZdt1+dZdt2)/2;
 case 'rk2',
  [dPdt,dZdt]=model(P,Z); P1=P+dt/2*dPdt; Z1=Z+dt/2*dZdt;
  [dPdt,dZdt]=model(P1,Z1);
  P=P+dt*dPdt; Z=Z+dt*dZdt;
 case 'rk4',
  [dPdt0,dZdt0]=model(P ,Z ); P1=P+dt/2*dPdt0; Z1=Z+dt/2*dZdt0;
  [dPdt1,dZdt1]=model(P1,Z1); P2=P+dt/2*dPdt1; Z2=Z+dt/2*dZdt1;
  [dPdt2,dZdt2]=model(P2,Z2); P3=P+dt  *dPdt2; Z3=Z+dt  *dZdt2;
  [dPdt3,dZdt3]=model(P3,Z3);
  P=P+dt*(dPdt0+2*dPdt1+2*dPdt2+dPdt3)/6;
  Z=Z+dt*(dZdt0+2*dZdt1+2*dZdt2+dZdt3)/6;
 case 'ab2',
  [dPdt1,dZdt1]=model(P   ,Z   );
  [dPdt2,dZdt2]=model(Pn_1,Zn_1);
  P=P+dt*(3/2*dPdt1-1/2*dPdt2);
  Z=Z+dt*(3/2*dZdt1-1/2*dZdt2);
 case 'ab2s',
  Pab=3/2*P-1/2*Pn_1;
  Zab=3/2*Z-1/2*Zn_1;
  [dPdt,dZdt]=model(Pab,Zab);
  P=P+dt*dPdt;
  Z=Z+dt*dZdt;
 case 'qab2',
  eps=0.1;
  [dPdt1,dZdt1]=model(P   ,Z   );
  [dPdt2,dZdt2]=model(Pn_1,Zn_1);
  P=P+dt*((3/2+eps)*dPdt1-(1/2+eps)*dPdt2);
  Z=Z+dt*((3/2+eps)*dZdt1-(1/2+eps)*dZdt2);
 case 'qab2s',
  eps=0.1;
  Pab=(3/2+eps)*P-(1/2+eps)*Pn_1;
  Zab=(3/2+eps)*Z-(1/2+eps)*Zn_1;
  [dPdt,dZdt]=model(Pab,Zab);
  P=P+dt*dPdt;
  Z=Z+dt*dZdt;
 case 'ab3',
  [dPdt1,dZdt1]=model(P   ,Z   );
  [dPdt2,dZdt2]=model(Pn_1,Zn_1);
  [dPdt3,dZdt3]=model(Pn_2,Zn_2);
  P=P+dt*(23*dPdt1-16*dPdt2+5*dPdt3)/12;
  Z=Z+dt*(23*dZdt1-16*dZdt2+5*dZdt3)/12;
 case 'ab3s',
  Pab=(23*P-16*Pn_1+5*Pn_2)/12;
  Zab=(23*Z-16*Zn_1+5*Zn_2)/12;
  [dPdt,dZdt]=model(Pab,Zab);
  P=P+dt*dPdt;
  Z=Z+dt*dZdt;
 otherwise,
  disp( sprintf('Error: method="%s" is not known',method) );
  error('Unrecongnized method!');
 end