%% set number of subintervals
Ni = [100 200 400 800 1600];

fprintf('N \t Its\n')
for i = 1:length(Ni)
    n = Ni(i);
    
    %% set up the matrix A
    e = ones(n, 1);
    x = linspace(0, 1, n+1)';
    h = x(2) - x(1);
    A = spdiags([e, -2*e, e], -1:1, n-1, n-1)/(h^2);
    
    %% Initial guess    
    u0 = 20*x.*(1-x);
    [u, its(i)] = nlnewton(A, u0(2:end-1), 0);
    
    %% pad boundary conditions in
    u = [0; u; 0];

    fprintf('%i \t %i\n', n, its(i))
    
    %% plotting both solutions at lowest grid resolution
    if i == 1
        hold on;
        u0 = zeros(size(x));
        [u, its(i), k] = nlnewton(A, u0(2:end-1), 1);
        u = [0; u; 0];
        plot(x, u, LineWidth=2);

        u0 = 20*x.*(1-x);
        [u, its(i), k] = nlnewton(A, u0(2:end-1), 1);
        u = [0; u; 0];
        plot(x, u, LineWidth=2);
    end
end

function [u, its] = nlnewton(A, u, pr)
    %% pr==1 enables printing and enables us to observe the quadratic 
    %% convergence of Newton's method
    its = 0;
    if pr
        fprintf('----------------------\n')
    end
    while 1
        %% Calculate Jacobian and function f
        J = A + diag(exp(u));
        f = A*u + exp(u);
        p = -J\f;
        its = its+1;
        if pr
            fprintf('it: %i    ||f||: %7.2e\n', its, norm(f));
        end
        
        %% Convergence criterion
        if norm(p) < 1e-8
            break
        end
        u = u + p;
    end
end
