function newton1D(x, f, dfdx, tol, maxit=100)
    
    for it in 1:maxit
        xnew = x - f(x)/dfdx(x)            
        if abs((xnew-x)/xnew) <= tol
            return xnew,it+1
        end
        x = xnew
        if it == maxit
            println("Warning, no convergence in $maxit iterations")
            return xnew,it+1
        end
    end
end;

import LinearAlgebra: norm

function multiNewt(f, Jac, x, tol, maxit=100)
    
    for it in 1:maxit
        fvec = f(x)
        J    = Jac(x)
        xnew = x - J\fvec            
        if norm(xnew-x)/norm(xnew) <= tol
            return xnew,it+1
        end
        x = xnew
        if it == maxit
            println("Warning, no convergence in $maxit iterations")
            return xnew,it+1
        end
    end
end;

fg(x) = [ x[1]^2 + x[2]^2 - 3,  
          2*x[2] - x[1]^2 ]

Jac(x) = [ 2*x[1]  2*x[2]; 
          -2*x[1]     2.0  ];

x0 = [10, 10]

x_solution,nit = multiNewt(fg, Jac, x0, 1.E-6)
println("solution: $x_solution")
println("f(x):     $fg(x_solution)")