clear all;

% Define function f and it's derivatives

a2 = 14;
a3 = 2;

f   =   @(x,a)([-x(1)+a*(1-x(1))*exp(x(2));
                -x(2)+a*a2*(1-x(1))*exp(x(2))-a3*x(2)]);
           
Dxf  =  @(x,a)([-1-a*exp(x(2)),   a*(1-x(1))*exp(x(2));
                -a*a2*exp(x(2)),  -1+a*a2*(1-x(1))*exp(x(2))-a3]);

Daf  =  @(x,a)([  (1-x(1))*exp(x(2));
                  a2*(1-x(1))*exp(x(2))]);


% Define derivative respect to u

fu  = @(u)( f(u(1:2),u(3)));
Dfu = @(u)( [Dxf(u(1:2),u(3)) Daf(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;0];
a = 0;
u0 = [ x0; a];

% initalize some stuff
%

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

% Continue with free parameters a1 and a2 from every Fold and Hopf
% use difference approximation for det( fx );
%

% define functions in case of free a1 and a2

ff   =  @(x,a)([-x(1)+a(1)*(1-x(1))*exp(x(2));
		-x(2)+a(1)*a(2)*(1-x(1))*exp(x(2))-a3*x(2)]);

Dxff =  @(x,a)([-1-a(1)*exp(x(2)),  a(1)*(1-x(1))*exp(x(2));
                -a(1)*a(2)*exp(x(2)),  -1+a(1)*a(2)*(1-x(1))*exp(x(2))-a3]);

Daff =  @(x,a)([   (1-x(1))*exp(x(2)) , 0 ;
		   a(2)*(1-x(1))*exp(x(2)) , a(1)*(1-x(1))*exp(x(2))]);

% define functions for Fold - bifurcations

gu = @(u)( [ff(u(1:2),u(3:4));
	    det(Dxff(u(1:2),u(3:4))) ]);

Dgu =@(u)([ Dxff(u(1:2),u(3:4)) , Daff(u(1:2),u(3:4));
	    Dapprox(@(x)(det(Dxff(x(1:2),x(3:4)))),u)]);

% Part of the Jacobian needed in Fold/Hopf monitoring
Dxgu =@(u)([ Dxff(u(1:2),u(3:4))]);

% continue from Fold points !

for i=1:size(Fold_point,2)
  
  % set a_2 = 14 !
  u0 = [ Fold_point(:,i) ; a2];
  
  % refine the initial point a littlebit 
  % u0 = gnewton(u0,gu,Dgu,@pifun);
 
  % do as previosly ( now in both directions )!
  
  [FU{i},FFP{i},FFI{i},FHP{i},FHP{i}] = path_follow(gu,Dgu,Dxgu,u0,0.05,5);
  [FUN{i},FFPN{i},FFIN{i},FHPN{i},FHPN{i}] = path_follow(gu,Dgu,Dxgu,u0,0.05,5,-1);
end



% define functions for Fold - bifurcations

gu = @(u)( [ ff(u(1:2),u(3:4));
	     det(bialt(Dxff(u(1:2),u(3:4)))) ]);

Dgu = @(u)([ Dxff(u(1:2),u(3:4))                Daff(u(1:2),u(3:4));
	     Dapprox(@(x)(det(bialt(Dxff(x(1:2),x(3:4))))),u) ]);


% continue from Hopf points !
for i=1:size(Hopf_point,2)
  
  % set a_2 = 14 !
  u0 = [ Hopf_point(:,i) ; a2];
  
  % refine the initial point a littlebit 
  %u0 = gnewton(u0,gu,Dgu,@pifun);
  
  % do as previosly ( now in both directions )!
  [HU{i},HFP{i},HFI{i},HHP{i},HHP{i}] = path_follow(gu,Dgu,Dxgu,u0,0.1,10);
  [HUN{i},HFPN{i},HFIN{i},HHPN{i},HHPN{i}] = path_follow(gu,Dgu,Dxgu,u0,0.1,10,-1);
  
end


% PLOT

norm_x =[];
a = [];
UU = U ;

for  i=1:1:size(UU,2)
  
  norm_x = [norm_x norm(UU(1:2,i))];
  a = [a UU(3,i)];
end

figure;plot(a,norm_x,'r-');

hold on;

for k=1:1:length(FU)
  UU = FU{k}
  norm_x = [];
  a = [];
  
  for i=1:1:size(UU,2)
    norm_x = [norm_x norm(UU(1:2,i))];
    a = [a UU(3,i)];
  end
  
  hold on;plot(a,norm_x,'b-');
  
end

for k=1:1:length(FUN)

  UU = FUN{k}
  
  norm_x = [];
  a=[];
  
  for i=1:1:size(UU,2)

    norm_x = [norm_x norm(UU(1:2,i))];
    a = [a UU(3,i)];
  end
  hold on;plot(a,norm_x,'g-');

end

for k=1:1:length(HU)
  UU = HU{k}
  norm_x = [];
  a = [];
  
  for i=1:1:size(UU,2)
    norm_x = [norm_x norm(UU(1:2,i))];
    a = [a UU(3,i)];
  end
  
  hold on;plot(a,norm_x,'b-');
  
end

for k=1:1:length(HUN)

  UU = HUN{k}
  
  norm_x = [];
  a=[];
  
  for i=1:1:size(UU,2)

    norm_x = [norm_x norm(UU(1:2,i))];
    a = [a UU(3,i)];
  end
  hold on;plot(a,norm_x,'g-');

end
