function [x,v,i] = newtcorr(x0, v0)
%
% Newton corrections, internal routine
%
global cds

x = x0;
v = v0;

R = [];
R(cds.ndim) = 1;

% WM: for-loops are faster then while-loops
for i = 1:cds.options.MaxCorrIters
  if i <= cds.options.MaxNewtonIters
    A = contjac(x);%x, A,pause
    B = [A; v'];
  end

  % repeat twice with same jacobian, calculating
  % the jacobian is usually a lot more expensive
  % than solving a system
  for t=1:2
    Q = [feval(cds.curve_func, x); 0];
    
    if cds.options.MoorePenrose
%      R = [A*v; 0];
      
      lastwarn('');
      D = B\[Q R'];
      %dx = B\Q; dv = B\R';
   
      if ~isempty(lastwarn)
        x = [];
        v = [];
        return;
      end
      
%      v = v - D(:,2);
      v = D(:,2);
      v = v/norm(v);
      
      dx = D(:,1);
    else
      dx = B\Q;
    end
    x = x - dx;
    
    % WM: funtol and vartol were switched
    if norm(dx) < cds.options.VarTolerance & norm(Q) < cds.options.FunTolerance
%      if ~cds.MoorePenrose 
        A = contjac(x);%,pause
        %    R = [A*v; 0];
        %    v = v - B\R;
        v = ([A;v']\R');
        v = v/norm(v);
%      end
      return;
    end
  end
end

x = [];
v = [];

%SD:Newton corrections pal/mp