function phaseplane

   global parameters

   a = 0.7;
   b = 0.8;
   c = 3.0;
   parameters = [a, b, c];

   u0 = [-1, -0.5]; 	% z = 0: small
%   u0 = [-1, -1]; 	% z = 0: large
   int_time = [0 200]; % Integration time.
   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 -1 1.75])
      hold on % "Hold" and wait for additional plotting information.  

      % NULL CLINES
      x= -3:0.05:3;
      y1= x-(1/3)*x.^3;
      y2= (a/b)+x/b;
      plot(x,y1,'m',x,y2,'m'); 
      z = 0.2;
      y11= x-(1/3)*x.^3+z;
      z = 0.5;
      y12= x-(1/3)*x.^3+z;
      plot(x,y11,'m--',x,y12,'m:');
      hold off % Stop holding and finish the plot.	
end

function dudt = f(t,u)

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

   if (50 < t) & (t < 100)
     z = 0;
     z = 0.2;
     %z = 0.5;
   else
     z = 0;
   end 

%   dudt(1) = c*(u(2)+u(1)-u(1)^3/3 + z);
%   dudt(2) = (1/c)*(a-u(1)-b*u(2));

   dudt(1) = c*(u(1)-u(1)^3/3 - u(2) + z);
   dudt(2) = (1/c)*(a+u(1)-b*u(2));


end
