clear all;

% Define function f and it's derivatives

f  =  @(x,c)([ (x(2)-x(2)^2-x(1)*( x(1)^2-x(2)^2+2*x(2)^3/3 + c));
	       (x(1)+(x(2)-x(2)^2)*(x(1)^2-x(2)^2+2*x(2)^3/3+c))]);

Dxf = @(x,c)([ -(x(1)^2-x(2)^2+2*x(2)^3/3+c)-2*x(1)^2, 1-2*x(2) - x(1)*(-2*x(2)+2*x(2)^2);
	      1 + (x(2) - x(2)^2)*2*x(1), (1 - 2*x(2))*(x(1)^2-x(2)^2+2*x(2)^3/3+c)+(x(2)-x(2)^2)*(-2*x(2)+2*x(2)^2)]);
	      
Dcf = @(x,c)( [ -x(1) ; x(2)-x(2)^2]);


% Define derivative respect to u

fu  = @(u)( f(u(1:2),u(3)));
Dfu = @(u)( [Dxf(u(1:2),u(3)) Dcf(u(1:2),u(3))] );
Dxfu = @(u)( Dxf(u(1:2),u(3)));


% use pseudo-arc length continuation to construct bifurcation
% diagram 
%
% Note that the needed trivial equilibrium point is c=0 , x=0
%
% Begin by computing the tangent of the curve i.e. vector from the
% nullspace of Duf
%

x0 = [0;1];
c = 0;
u0 = [ x0; c];


delta = 0.2;

[U, Fold_point, Fold_indicator, Hopf_point, Hopf_indicator] = path_follow(fu,Dfu,Dxfu,u0,0.2,10)