function phaseplane

   global eps beta	
   
   eps = 0.005;
   beta = 5.13;
   u0 = [1.,  2.0,0]; 	

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

   figure(1) 
    subplot(3,1,1)
    plot(t,u(:,1),'b');
    xlabel('time')
    subplot(3,1,2)
    plot(t,u(:,2),'k');
    subplot(3,1,3)
    plot(t,u(:,3),'g');


    figure(2)
      plot(u(2000:end,3),u(2000:end,1),'b');
      xlabel('u(3)=c')
      ylabel('u(2)=u')   
      axis([-0.5 1.25 -2.5 2.5]) 
      hold on 
      % Bifurcation equation
      x= -3:0.05:3;
      y1= 1 - x.^3-2*x.^2;
      y2= 4*x+beta;
      plot(y1,x,'m',y2,x,'m'); 
      hold off 
end

function dudt = f(t,u)

   dudt = zeros(3,1); % Create a column vector.

   global eps beta	
   
   dudt(1) = u(2) - u(1)^3 + 3*u(1)^2 - u(3);
   dudt(2) = 1 - 5*u(1)^2-u(2);
   dudt(3) = eps*(u(1)+(beta - u(3))/4);

end
