function [x,fv]=NewtonLSQ(x,xi,eta)
% NEWTONLSQ Curve fiting with Newton's method 
%   [x,fv]=NewtonLSQ(x,xi,eta) fits the function f(t)=x1+x2*exp(-x3*t)
%   to given points (xi,eta). x is the initial guess and is overwritten
%   with fitted parameters x1,x2 and x3. fv contains norm(f) for
%   each iteration

h=x; fv=[];
while norm(h)>100*eps*norm(x)
  ee=exp(-x(3)*xi);
  tee=xi.*ee;
  J=[ones(size(xi)),ee,-x(2)*tee ];          % Jacobian
  f=x(1)+x(2)*ee-eta;                        % residual
  s1=f'*tee; s2=x(2)*f'*(xi.*tee);
  A=J'*J+[0 0 0; 0 0 -s1; 0 -s1 s2];         % Hessian of Phi
  h=-A\(J'*f);
  x=x+h;
  xh=[x,h]
  res=norm(f)
  fv=[fv norm(f)];
end