function [x1, n] = newtonsys (f, J, x0, eps, nmax)
    dx = J(x0)\-f(x0); % first iteration
    x1 = x0 + dx;       % first iteration
    n = 1;
    while norm(x1-x0)>eps && n<=nmax
        x0 = x1;
        dx = J(x0)\-f(x0);
        x1 = x0 + dx;
        n = n + 1;
    end;
end