function predprey

   global parameters % Share with subroutines/functions

   % predprey parameters
   parameters = [0.4, -0.3, -0.01, 0.005]; % [r1,r2,a1,a2]

   int_time = [0 200]; % Integration time.
   u0 = [60, 60]; % Initial conditions
   options = odeset('AbsTol',1e-8,'RelTol',1e-8); % Control accuracy of solver
   [t,u] = ode45(@f,int_time,u0,options);

   figure(1) % Solutions as a function of time.
   subplot(2,1,1) % 2x1 array of plots, 1st plot.
      plot(t,u(:,1),'b',t,u(:,2),'k');
      xlabel('time')
      ylabel('u(1) blue, u(2) black')

   subplot(2,1,2) % 2x1 array of plots, 2nd plot.
      plot(u(:,1),u(:,2));
      xlabel('u(1)')
      ylabel('u(2)')

      % VECTOR FIELD: Comment out if not needed.	
      hold on % "Hold" and wait for additional plotting information.  
      x = linspace(0, 120, 20); % Values of u1 to plot vectors.
      y = linspace(0,  80, 20); % Values of u2 to plot vectors.
      [u1grid,u2grid] = meshgrid(x,y); % Create 2 matrices that hold the values of u1 and u2 at each point. 
      Du1 = zeros(size(u1grid)); % Matrix that will hold the derivative of u1.
      Du2 = zeros(size(u2grid)); % Matrix that will hold the derivative of u2.
      t = 0;
      for i = 1:numel(u1grid)
         uderivatives = f(t,[u1grid(i);u2grid(i)]);
         Du1(i) = uderivatives(1);
         Du2(i) = uderivatives(2);
      end
      quiver(u1grid,u2grid,Du1,Du2,'r'); % Plot the vectors [Du1,Du2] at the grid points [u1grid, u2grid].
      axis tight; % Set the axis range to the data range.
      hold off % Stop holding and finish the plot.

end

function dudt = f(t,u)

   global parameters
   dudt = zeros(2,1); % Create a column vector.

   r1 = parameters(1);
   r2 = parameters(2);
   a1 = parameters(3);
   a2 = parameters(4);
	
   dudt(1) = r1*u(1) + a1*u(1)*u(2);
   dudt(2) = r2*u(2) + a2*u(1)*u(2);

end
