function phaseplane

   global a

   a = 0.99645152
   u0 = [1, 1]; 
   int_time = [0 100]; 
   options = odeset('AbsTol',1e-8,'RelTol',1e-8); 
   [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([-8 8 -4 4])

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

% NULL CLINES
      
      y1= a;
      y2= -8:0.05:8;
      x = -2*(y2.^3-3*y2);
      plot([-8,8],[a,a],'m',x,y2,'m'); 
      hold off % Stop holding and finish the plot.	

      hold off % Stop holding and finish the plot.	
end

function dudt = f(t,u)

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

end
