function x = newton_solver(f,m,x0,dfdx)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                              %
% function x = newton_solver(f,m,x0,dfdx)  %
%                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                              %
% Solves the equation f(x) = m using x0 as     %
% start guess where f is a string with the name%
% of the file in which f is defined.           %
% The jacobian can be specified                %
% in the file 'dfdx'. If omitted, a numerical  %
% approximation is used                        %
%                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


tol=0.01;
x=x0;
fx=feval(f,x0);
err=1;
errmax=100; 

if nargin == 4
% Use exact jacobian
while (err>tol)
    fx = feval(f,x); 
    J = feval(dfdx,x);
    dx = -J\(fx-m);
    x = x+dx;
    fx = feval(f,x);
    err=norm(fx-m,'inf');
    if(err>errmax)
        error('Divergent Newton solver (file: newton_solver.m)')
    end
end
else
% Use approximate jacobian
while (err>tol)
    fx = feval(f,x);
    J = appr_dfdx(f,x);
    dx = -J\(fx-m);
    while (norm(feval(f,x+dx)-m)>norm(fx-m))      % added 000211
        dx=dx./2
    end						  % end added
    x = x+dx;
    fx = feval(f,x);
    err=norm(fx-m,'inf');
    if(err>errmax)
        error('Divergent Newton solver (file: newton_solver.m)')
    end
end
end