%% Solving an ODE (NPZD) 
% This script shows you how to solve a system of ODEs (the NPZD model in
% this case). It uses the built in Matlab function ode45 which uses an
% explicit Runga-Kutta (4,5) method. It has an adaptive time step so time
% is not equally spaced. Please note that this script defines functions at
% the end, which is only supported by MATLAB 2016b or later.
% Marek: it would be easy to make a back compatible version of this by
% using inline functions.
%% Initialization
% We begin by defining all the parameters in the NPZD model. Since we want
% to use them inside of the function that defines the ODE we define them to
% be "global variables" (this means that they can be used inside of a
% function which normally has a separate scope).
global kn mu mp kp k Im alpha mz beta
mu = 2;
kn = 1;
mp = 0.1;
mz = 0.2;
alpha = 0.7;
beta = 0.3;
kp = 0.7;
Im = 1.5;
k = 0.1;
%%
% We define the ICs in the order: N; P; Z; D.
V0 = [1.4; 0.2; 0.2; 0.2];
%% Running the solver
tic
% To use the ode45 function you need to provide: the function that defines
% the DE, a vector indicating the initial time and final time, and the ICs
% of your variable. Please see the npzd function at the end for more about
% the DE.
[t, V] = ode45(@func_npzd,[0 60],V0);
toc
plot(t,V)
%%
% You can try changing the various parameters and see how the system
% changes. Some interesting parameters to change are the ICs, k and mu.
%% Notes
% ode45 is the standard ODE solver and will work for most cases. If you
% have a stiff ODE system, this method will not work and you will have to
% use another method such as ode23s. To see more about the different built
% in Matlab ODE solvers check out: https://www.mathworks.com/help/matlab/ordinary-differential-equations.html

%% NPZD ODE
% NPZD is the differential equation for modeling Nutrient, Plankton,
% Zooplankton and Detritus.
function dvdt = func_npzd(t,V)
%%
% All variables are contained in the vector V and their derivatives are
% returned in dvdt. Requires that the parameters have been defined
% globally.
%% Explicitly declare names
% We define each variable explicitly to make it easier to write the
% equations
N=V(1);
P=V(2);
Z=V(3);
D=V(4);
%% 
% Call the globally defined variables to have them defined in our scope.
global kn mu mp kp k Im mz alpha beta
%% Equations
% The set of DEs in the form dY/dt = f(Y);
dPdt = (N/(kn + N))*mu*P - mp*P-(P/(kp+P))*Im*Z;
dZdt = alpha*(P/(kp+P))*Im*Z-mz*Z;
dDdt = mp*P + (1-alpha-beta)*(P/(kp+P))*Im*Z-k*D;
dNdt = -dPdt-dZdt-dDdt;
%% Return
% Pass the derivatives back
dvdt = [dNdt; dPdt; dZdt; dDdt];
end
