%
% chap2page39.m
%
% second order finite difference method for the nonlinear ODE bvp
%   u''(x)+sin(u(x))=0, ax<x<bx
%   u(ax)=alpha, u(bx)=beta
% Using central differences on a uniform grid.
%
% Solve ODE (2.77) with Newton's method (2.80-2.82)
% Based on  http://www.amath.washington.edu/~rjl/fdmbook/chapter2  (2007)
set(0, 'defaultaxesfontsize',16,'defaultaxeslinewidth',1,...
    'defaultlinelinewidth',1,'defaultpatchlinewidth',0.7,...
    'defaulttextfontsize',16);

%setup the ODE
ax = 0;
bx = 2*pi;
alpha = 0.7;
beta = 0.7;
%utrue = @(x) exp(x);  % true soln unknown

%construct f(x) based on the ODE
f = @(x) zeros(size(x));  % right hand side function


% Solve the problem for ntest different grid sizes to test convergence:
m1vals = [100];
ntest = length(m1vals);
hvals = zeros(ntest,1);  % to hold h values
E = zeros(ntest,1);      % to hold errors

figure; hold on
for jtest=1:ntest
    m1 = m1vals(jtest);
    m2 = m1 + 1;
    m = m1 - 1;
    h = (bx-ax)/m1;
    hvals(jtest) = h;
    x = linspace(ax,bx,m2)';   % m+2 points counting the boundaries
    
    kmax=5; %max iteration
    err=1; tol=1e-10;
    %u0=0.7*cos(x(2:end-1))+0.5*sin(x(2:end-1));%initial guess, Fig 2.4(a)
    %u0=0.7*ones(size(x(2:end-1)));%initial guess, Fig 2.4(b)
    u0=0.7+sin(x(2:end-1)/2);%initial guess, Fig 2.5
    %Start Newton iteration as inner loop
    k=0;
    while (err > tol && k <kmax)
      
        plot(x,[alpha;u0;beta])
        % set up Jacobian matrix J (using sparse matrix storage):
        % construct system (2.82) by first compute 3 main diagonals
        %
        e=ones(size(x(2:end-1)));
        J = 1/h^2 * spdiags([e -2+h^2*cos(u0) e], [-1 0 1], m, m);
        %construct the nonlinear (2.78) system
        A = 1/h^2 * spdiags([e -2*e e], [-1 0 1], m, m);
        % Right hand side: see (2.78)
        F = f(x(2:end-1)); %only interial points
        F(1) = F(1)-alpha/h^2;
        F(m) = F(m)-beta/h^2;
        Gu=A*u0+sin(u0)-F; %(2.78)
        % solve linear system:
        deltaU = J\(-Gu);
        unew=u0+deltaU;
        u0=unew; %update solution
        % compute the difference between two iterations
        k=k+1;
        err=max(abs(deltaU))   
        
    end
    
    
    % estimate errors from the difference between two iterations
    %E(jtest) = max(abs(deltaU));
end

%error_table(hvals, E);   % print tables of errors and ratios
%error_loglog(hvals, E);  % produce log-log plot of errors and least squares fit

%plot the solutions and the error curves


