% mnewton.m: Newton's method for a system of nonlinear equations
% Usage: Define subroutines f and df, call as mnewton(@f, @df, m, n, x0, eps)
% Here "eps" is tolerance for root, x0 is initial guess. Returns a vector 
% "r", the solution. Here m is number of equations, n number of vars.
function r = mnewton(f, df, m, n, x0, eps)

[p,q] = size(x0);
if(p==1)
    r = x0';
else
    r=x0;
end

xold = r+eps; %prevent termination on first iteration

while(norm(r-xold)>eps)
    
    xold = r;
    J = df(r,m,n);
    fx = f(r,m);
    r = r-J\fx
    
end

