clear all

n=99
dx=1/(n+1)
newtontol=10^(-10)
gamma=0.9
gamma_max=0.999

for i=1:n
    x(i)=i/(n+1);
    f(i)=2-sin(x(i)*(1-x(i)));
end
f=f';
x=x';

A = 2*eye(n) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1);
A = A/dx/dx;

rhs = f;
newtoniter=0;
gmresiter=0;
newton_res=norm(rhs)
convtest = newton_res;
newtontol = newtontol*newton_res
uvec=zeros(n,1);

while (newton_res>newtontol)
    if (newtoniter==0)
        newton_res_old = newton_res;
	    linrtol  = gamma_max
    else
        gammaA = gamma*newton_res^2/newton_res_old^2;
        if (gamma*linrtol^2 < 0.1)
            gammaC = min(gamma_max,gammaA);
        else
            gammaC = min(gamma_max, max(gammaA,gamma*linrtol^2));
        end
        linrtol = min(gamma_max,max(gammaC,0.5*newtontol/newton_res))
        newton_res_old = newton_res;
    end
    M = A-diag(cos(uvec));
    linrtol=10^(-10);
    [deltau, flag, relres, iter] = gmres(M,rhs,[],linrtol, n);
    iter
    uvec = uvec +deltau;
    rhs = -A*uvec+sin(uvec)+f;
    newton_res=norm(rhs)
    newtoniter = newtoniter+1
    convtest = [convtest; newton_res];
    gmresiter = gmresiter+iter;
end

plot(x, uvec)

newtoniter
gmresiter
