function phaseplane


   int_time = [0 20]; % Integration time.
%>>>> Edit to get 3rd example.
   u0 = [1, 1]; % ICs for 1st and 2nd example
   % u0 = [-10,-10]; % ICs for 3rd example
   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)')
      axis([-2 2 -2 2])

      % VECTOR FIELD: Comment out if not needed.	
      hold on % "Hold" and wait for additional plotting information.  

%>>>> Edit to get 3rd example.
      limits = [-2,2]; %Limits for 1st and 2nd example.
      %limits = [-0.1,4]; %Limits for 3rd example

      x = linspace(limits(1), limits(2), 20); % Values of u1 to plot vectors.
      y = linspace(limits(1), limits(2), 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.

      % NULL CLINES
      % Plot lines by specifying the range of u1 and the range of u2.
      % plot( [u1low, u1high],[u2low,u2high],'magenta');
      plot([-4,4],[0,0],'m'); % u2 = 0.
      plot([0,0],[-4,4],'m'); % u1 = 0.
%>>>> Edit to get 3rd example.
      %plot([-4,4],[1,1],'g'); % u2 = 1. Additional NC for 3rd example.
      %plot([1,1],[-4,4],'g'); % u1 = 1. Additional NC for 3rd example.  

      % Given an equation for a nullcline u2 = g(u1).
      % Generate an array of u1 values.
      %   u1vals = u_start:step:u_stop;
      % Substitute them into the function g. 
      %    u2vals = g(u1vals);
      % Put a . before multiplication, division or power operations (other than constants) 
      % so that matlab knows to operate element wise. For example, if u2 = u1^2/(1+2*u1). 
      %    u2vals = u1vals.^2./(1 + 2*u1vals);     

      hold off % Stop holding and finish the plot.	
end

function dudt = f(t,u)

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

   %dudt(1) = -u(1);
   %dudt(2) = -3*u(2);

%>>>> Edit to get 3rd example.
   %dudt(1) = -u(1) + u(1)*u(2);
   %dudt(2) = +u(2) - u(1)*u(2);
   
   %dudt(1) = u(1)^2 -1;
   %dudt(2) = u(1)*u(2);

%>>>> Limit cycles
   w=1;
   s=1;   
   dudt(1) = -w*u(2) + u(1)*(s-(u(1)^2+u(2)^2));
   dudt(2) =  w*u(1) + u(2)*(s-(u(1)^2+u(2)^2));

%>>>> Modified van der pol oscillator
%        d = 1;
%        dudt(1) = u(2);
%        dudt(2) = -d*(u(2)^3 - u(2)) - u(1);




end
