% This program solves a modified version of the Lotka-Volterra predator-prey system of differential equations. The nondimensional model has the form
% dn/dt = n(1-n) - anp/(n + d)
% dp/dt = bp(1-p/n)
% where n = prey, p = predator, and a,b,d are positive parameters.
% We are going to set up a solution matrix, u, whose columns will be the solutions of n and p at various times, t.
% Because we have two equations, the matrix u will have two columns and many rows, depending on how many time steps there are.
%The more equations you have, the more columns the solution matrix will have.
function predprey % Define a function called predprey, note that this function has the same name as
% the .m file.
tspan = [0; 500]; %This command defines the time interval of interest.
u0 = [.4 .35]; %This command tells MATLAB the initial condition. Note that u0 is a (1 x 2) vector
%whose first entry is n(0) and second entry is p(0). If you are solving a system of more than
%two equations, simply place the additional initial conditions into the vector u0.
[t,u] = ode15s(@f,tspan,u0); % This command tells MATLAB to use the differential equation solver called ode15s to
% numerically compute the solution of the system of differential equations defined by the
% function f (below), for time interval tspan, and initial condition u0. The right-hand
% side of the equation tells MATLAB to store this output in an array containing
% a vector t (for the time points) and a matrix u (for the population densities).
%This command will not change even if you are solving a larger system of equations.
for i = 1:length(t) %This for-loop extracts the solutions for each variable, n and p, from the matrix u.
n(i,1) = u(i,1); %Specifically, the first column of u is the prey population, n and the second column of u is p.
p(i,1) = u(i,2); %If you are solving a larger system of equations, simply continue to name the additional columns
end %of the u matrix. Eg: m(i,1) = u(i,3) and so on.
figure;
%plot(n,p); %This tells MATLAB to plot the solution in the phase plane.
plot(t,n,t,p) %This tells MATLAB to plot the solutions as a function of time.
% --------------------------------------------------------------------------
function dudt = f(t,u) %This commands defines the function du/dt = f(t,u). Remember that u is a matrix that contains both n and p.
a = 1; %Now define the parameters.
b = .05;
d = 0.2;
dudt = [u(1)*(1-u(1)) - a*u(1)*u(2)/(d+u(1)); b*u(2)*(1-u(2)/u(1))]; %This command inputs the left-hand side of the system of differential equations.
%Note that each equation is separated by a semicolon and here the n-variable is
%denoted u(1) and p is u(2). If you had more equations, you would input them here %and separate them by a semicolon. The third variable would be u(3) and so on.