function gn
clear all; 
format compact; 
format long;
%Set start value 
x0=...                       % case (i)
% x0=...                     % case (ii)
%Set Tolerance values 
rtol=1e-8; epsilon=1e-12;
%With minimisation
withmini=1;
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,@funnew)
%Without minimisation
withmini=0;
[iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,@funnew)


function [iter,x,resid] = gaussnewton(x0,rtol,epsilon,withmini,f)
% Gauss-Newton for nonlinear overdetermined system, 
% Numerische Mathematik, D-MAVT FS 13, K.Nipp, Serie 10, exercise 3
%
% in:   x0       starting values for iteration
%       rtol     stop if ||x-xold|| <= rtol*||x||
%       epsilon  stop minimisation if F(y) < F(x)+epsilon
%       withmini ONLY do the minimisation if set to 1
%       f        function containing the evaluation of f(x)-z
%                and its Jacobian. One output argument gives
%                just the residual
% out:  iter     number of iterations needed till convergence
%       x        solution in least square sense
%       resid    2-norm of resulting residual (for final x)

x=x0; 
xold=x0+1; 
iter = 0;
while (norm(x-xold)>rtol*norm(x) && iter<1000)
  iter = iter + 1; xold = x;
  % solve C*xi=b (C=C(x), b=-d=-(f(x)-z)) in least square sense
  [rx,C]=feval(f,x);        % residual at x
  xi = C \ -rx;  % C*xi=-rx
  if withmini==1  % minimisation [Schwarz sec. 7.4.2]
    tt = 1; Fx = rx'*rx;
    while (1)
      y = x + tt*xi; ry = feval(f,y); Fy = ry'*ry; 
      if (Fy < Fx+epsilon) 
          x = y; 
          break; % accept y as next x
      else tt = tt/2; % reduce tt and try again 
      end
    end
  else x = x + xi; % 'normal' newton (no minimisation, i.e. tt=1)
  end
end
resid=norm(feval(f,x));
if (iter==1000 || ~(resid<=10)) 
    error('no convergence'); 
end
%-------------------------------------------------------------------------
function [F,J]=funnew(x)
z=[3.85 2.95 2.63 2.33 2.24 2.05 1.82 1.8 1.75]';
t=[0 0.5 1 1.5 2 3 5 8 10]';
zt = %!! z(t)
F=zt-z;
if nargout > 1
J = %!! Jacobi-Matrix J(x) von F
end
