function [x, k] = newtonsystem(f, jac, x0)
% NEWTONSYSTEM Newton's method for a system of nonlinear equations.
% Input:
%   f        function that computes residual (n-vector)
%   jac      function that computes Jacobian matrix
%   x0       initial root approximation (n-vector)
%
% Output       
%   x        array of approximations (one per column, last is best)
%   k        number of iterations

    % Operating parameters.
    funtol = 1000*eps;  
    xtol = 1000*eps;  
    maxiter = 40;

    x = x0(:);  
    y = f(x0);
    J = jac(x0);
    dx = Inf;
    k = 1;

    while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
        dx = J\y;   % Newton step
        x(:,k+1) = x(:,k) - dx;
    
        k = k+1;
        y = f(x(:,k));
        J = jac(x(:,k));
    end
end
