function nonlin_gmres_demo()
fprintf(1,'\n*** Small nonlinear system:\n');
n=15;
u=nonlin_gmres(zeros(n,1),20);
h=1/(n-1);
n=150;
fprintf(1,'\n*** Bigger nonlinear system:\n');
uu=nonlin_gmres(zeros(n,1),20);
hh=1/(n-1);
%
% plot solutions of the differential equation
%
plot(0:h:1,u,0:hh:1,uu);


function u=nonlin_gmres(x0,itmax)
global x
x=x0;
for iter=1:itmax
    f = resid(x);
    dx = gmres(@afun,f,10,1e-3,500);
    x = x - dx;
    fprintf(1,' Newton iteration #%d: |dx|=%.2e\n',iter,norm(dx));
    if norm(dx) < 1e-4
        u = x;
        return;
    end
end
u=x;
fprintf(1,'Itmax exceeded...\n');


function y=afun(v)
global x
ep = 1e-3;
y = (resid(x+ep*v)-resid(x))/ep;

%
%  Nonlinear system R(X)=0 corresponding discretization of
%  nonlinear ordinary differential equation:

%  -u''(x)=exp(u(x))   0<x<1
%   u(0)=u(1)=0
%
function r = resid(x)
n=length(x);
r=zeros(n,1);
h=1.0/(n-1);
r(1)= x(1);
r(n)= x(n);
for i=2:n-1
   r(i) = 2*x(i)-x(i+1)-x(i-1)-h*h*exp(x(i));  
end
