function [x, eps, info] = newrap(f, x, tol)
%% ... Newton Raphson method for solving a system of
%% ... nonlinear equations. Jacobi matrix is evaluated
%% ... using central differences
%% ... http://www.pcs.cnu.edu/~bbradie/msystems.html
%% check sizes
[neq,trs] = size(x);
J         = zeros(neq,neq);
w         = zeros(neq,1);
fx        = zeros(neq,1);
delta     = zeros(neq,1);
eps       = zeros(neq,1);
niter = 0;
h     = 1e-6;
eps   = f(x);
while (abs(eps) > tol & niter < 100)
    niter = niter + 1;
    
    %% ... evaluate jacobian - forward diff.
    for i = 1 : neq
        %% ... f returns a row vector fi(x1,...,xN)
        %% ... thus the deriv. ---> df(1)/dx(i), ... , df(N)/dx(i) 
        %% ... ! Improve this! - Richardson extrapol. ?
        a = x;
        a(i) = x(i)+h;
        J(:,i) = (f(a)-f(x))/(h);
    end
    
    %% solve J*delta = -f(x)
    delta = J\(-f(x));
    x     = x + delta;
    eps   = f(x);
    %J
    %x
    %epsilon
end % while
info = niter;
return;
end %%function