% Script for Newton's method on nonlinear pendulum BVP
%          d2u/dt2 + sin(u) = 0,  u(0) = a, u(T) = b
% on a grid with t(j) = jT/N, j=1,...,N-1. 
% For this example, we use a grid spacing h = T/40, a = 0.7, b = 0.7, T = 2*pi
% and an initial guess u(t) = a for all t.

%  Define parameters

  a = 0.7;
  b = 0.7;
  T = 2*pi;

  N = 40;  % Number of gridpoints
  h = T/N;
  m = N-1;
  tol = 1e-6;  % If max norm of diff between iterates < tol, assume convergence
  maxiter = 20; % Maximum number of Newton iterations

%  Initial guess

  u = a*ones(m,1);

%  Newton iteration.  Allow a maximum of maxiter iterations to prevent infinite
%  loop.

  clf
  subplot(2,2,1)
  
  for iter = 1:maxiter
    
    % Plot solution before iteration
    plot((0:N)*h,[a;u;b]','g--')
    hold on

    G = diff([a;u;b],2)/h^2 + sin(u);   % Objective function vector
    e = ones(m,1);
    J = spdiags([e/h^2 -2*e/h^2+cos(u) e/h^2], -1:1, m, m);%J_ij = dG(i)/du(j) 
    du = -J\G;
    u = u + du;
    dumax = norm(du,inf)    % Display norm of du at each iteration
    if (dumax < tol)        % Exit loop when norm of du < tol.
      break
    end
  end

%  Now plot converged solution vs. t

  plot((0:N)*h,[a;u;b]')
  xlabel('t')
  ylabel('u')
  axis([0 2*pi -1 1])
  title('Nonlinear Pendulum BVP')
  hold off
