clear all;

% Define all needed functions here !

% helper function - G(x)
G    = @(x)( [ 10*(x(2)-x(1)^2);
              1-x(1);
	      3*sqrt(10)*(x(4)-x(3)^2);
	      1-x(3);
	      sqrt(10)*(x(2)+x(4)-2)
	      (x(2)-x(4))/sqrt(10)]);

% Jacobian of G
DXG  = @(x)( [ -20*x(1) 10 0 0 ;
	      -1       0  0 0;
	      0        0 3*sqrt(10)*[-2*x(3), 1];
	      0        0 -1 0;
	      0        sqrt(10)*[1 0 1];
	      0       [1 0 -1]/sqrt(10)]);

% Derivatives of the DG
DXX1G = [ -20 0 0 0; zeros(5,4)]; 
DXX2G = zeros(6,4);
DXX3G = [ zeros(2,4); 0 0 -6*sqrt(10) 0; zeros(3,4)]; 
DXX4G = zeros(6,4);

% Function g(x) and its Jacobian
g = @(x)(DXG(x)'*G(x));
Dg =@(x)([DXX1G'*G(x) ,DXX2G'*G(x), DXX3G'*G(x),  DXX4G'*G(x)]+ DXG(x)'*DXG(x) );

% Initial Guess
x0 = [-3 -1 -3 -1]';

% function h(x,a) goes into the predcorr
h   = @(x,a)( g(x) +(a-1)*(g(x0)) );
DXh = @(x,a)( Dg(x) );
DAh = @(x,a)( g(x0) );

% solution at a = 0
x = x0
a = 0;

% Max iteration and initial Delta
maxIter = 1000;
Del = 0.2;

% Save iterates into these variables
X = [x];
A = [0];

% Select the first corrector-direction as zero. The predcorr -
% corrects the direction.
tan = zeros(4,1);

for i=1:1:maxIter
  
  % use predcorr2 ( More efficent )
  [x,a,Del,tan] = predcorr2(x,a,Del,tan,h,DXh,DAh);

  % save iterates
  X = [X x];
  A = [A a];
  
  % If needed, modify the step, so that a = 1 will be hit
  if( (a+Del) > 1)
    disp('Del modified');
    Del = 1-a;
  end
  
  % Check if completed
  if( a == 1)
    break;
  end
  
end


% Plot error etc.. figures

figure;
for i=1:1:size(X,2)
  error(i) = norm(g(X(:,i)));
end

semilogy(1:size(X,2),error,'r*');
