function [x0,y0,err,His] = newton_mv(f,jac,x0,delta,epsilon,max1) 
%---------------------------------------------------------------------------
%Newton's method.

x0=x0(:); % column vector
His = x0';
y0 = feval(f,x0);
for k=1:max1,
  J = feval(jac,x0);
  [L,U,P]=lu(J); 
  if min(abs(diag(U))) == 0,
    dx = 0;
  else
    dx = -U\(L\P*y0);  % fazer a subst. e retros.
  end
%  dx=-J\y0;
  x1 = x0 + dx;
  y1 = feval(f,x1);
  err = norm(dx);
  relerr = err/(norm(x0)+eps);
  x0 = x1;
  y0 = y1;
  His = [His;x1'];
  if (err<delta)&(relerr<delta)&(norm(y1)<epsilon), break, end
end
