% 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.	
												  	
