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 < 1000)
    niter = niter + 1;
    %% ... evaluate jacobian - forward diff.
    [J,err] = jacobianest(f,x);
    %%[J,err] = derivest(f,x,'MaxStep',1000,'vectorized','yes')
    %% solve J*delta = -f(x)
    delta = J\(-f(x));
    x     = x + delta;
    eps   = f(x);
end % while
info = niter;
return;
end %%function