
% define needed functions
Uf = @(u)( [ -u(1)+u(3)*(1-u(1))*exp(u(2)) 
	      ; -3*u(2)+14*u(3)*(1-u(1))*exp(u(2))]);

DUf = @(u)( [ -1-u(3)*exp(u(2)), u(3)*(1-u(1))*exp(u(2)),(1-u(1))*exp(u(2)) ; 
		-14*u(3)*exp(u(2)) , -3+14*u(3)*(1-u(1))*exp(u(2)), 14*(1-u(1))*exp(u(2))]);


% the f - function staisfies f(0,0) = 0.
% this is the point, where we start the pseudo-arc-continuation.
x = [0;0];
a = 0;

u = [x;a];

% save the point
U = [u];

% constants !
maxIter = 1000;
Delta = 0.01;
 
% for the first step select v0 from the N(fx(x)).  
% use the following :  Q*R*x = 0 => R x = 0 , if x != 0
% note, that eig - normalizes all eigenvetors => |v| = 1

[Q,R] = qr( DUf(u));
[X,D] = eig( [R ; zeros(1,3)]);

I = find(abs(diag(D))<1e-10);
v = X(:,I(1));

while u(3) < 1;

  % take a Delta step  
  uold = u;
  u= u + Delta*v;
  
  % iterate with Newton and QR-dec.( use non-stationary Newton )
  u = gnewton(u, Uf, DUf);
 
  % construct new direction 
  v = (u-uold)/norm(u-uold);

  if( u(3) > 1+1e-8)
    u = uold;
    Delta = Delta/2;
  end
  
  U = [U u];
  
end

for i=1:length(U)
  nx = norm(U(1:2,i))
end;

figure;
plot(U(3,:),nx);

figure;
plot(U(3,:),U(1,:));

figure;
plot(U(3,:),U(2,:));

figure;
plot3(U(1,:),U(2,:),U(3,:));