function x=nsolve(a,b,x0,varargin)
% solve nonlinear system a(x)=b using Krylov iteration with
% Jacobian-vector multiplication by finite differences
% extra arguments are passed to a
% arguments
% a         the function
% b         the rhs
% x0        initial approximation
% varargin  optional extra arguments for a
% output
% x         approximate solution

newt_tol   = 1e-5;
newt_maxit = 5;
h = 1e-1;
krylov_tol   = sqrt(newt_tol);
krylov_maxit = min(10,length(b));
x = x0;

for it=1:newt_maxit
    % one Newton step
    % solve the linear system a(x)+a'(x)*dx=b
    ax = feval(a,x,varargin{:});
    % get the rhs of the linear system
    r = b-ax;
    err=norm(r);
%    fprintf('Newton iteration %i starting error %g\n',it,err)
    if err<newt_tol, break, end
    % solve the system by GMRES iterations
    switch version('-release')
        case '14'
        gm=@gmres_r14_jmfix;
        otherwise
        gm=@gmres;
    end
    [dx,tmp]=feval(gm,@amult,r,[],krylov_tol,krylov_maxit,[],[],[],a,x,ax,h,varargin{:});
    relres = norm(amult(dx,a,x,ax,h,varargin{:})-r,inf)/norm(r,inf);
%    fprintf('Jacobian system relative residual %g\n',relres)
    % newton update
    x = x+dx;
end   
    
function ad = amult(d,a,x,ax,h,varargin)
% matrix-vector multiplication by jacobian of a from finite differences around x
ad=(feval(a,x+d*h,varargin{:})-ax)/h;