% Pendulum BVP
%
% Dirichlet boundary conditions
%
global T alpha beta
%
T=2*pi;
alpha=0.7;
beta=0.7;
%
%hold off
mm=2^3;
hh=[];
tol=1e-8;
%%%%%%%%%%%%%%%%%%%
% computation of a reference solution
% we will use it in place of the exact solution
% we use the difference method with a very
% fine grid, i.e. h=1/(2^15), to compute the reference solution
% we use a smaller tolerance, tol/1000, in the Newton iteration, in the 
% computation of the reference solution
%%%%%%%%%%%%
changeinval=0;
display('choose the initial guess for the Newton iteration')
display('different initial guesses lead to different solutions of G(theta)=0')
display('"1" is 0.7+sin(x/2), "2" is 0.7cos(x), ">2" is 0.7cos(x)+0.5sin(x)')
changeinval=input('1,2,3,?')

mm=mm*2^12;
h=T/(mm);
m=mm-1;
xx=[h:h:T-h]';
theta=zeros(m,1);
if (changeinval == 1 )
   theta0=0.7*ones(m,1)+sin(xx/2);%
   else if (changeinval == 2)
   theta0=0.7*cos(xx);
        else
        theta0=0.7*cos(xx)+0.5*sin(xx);
   end
end   
%
theta=theta0;
maxNit=100;
%
% Newton iteration to solve G(theta)=0
%
    for k=1:maxNit
       thetaold=theta;
       J=JacG(theta);
       g=G(theta);
       theta=theta-J\g;
       % difference between two successive iterations as stopping criterion
       nr=norm(theta-thetaold)/norm(theta0)
       if (nr <tol/1000 )
           display('REFERENCE SOLUTION')
           display('Number of Newton iterations to satisfy the convergence criterion')
           k
           break
       end    
    end
    k
thetaex=theta; 
xxx=xx;
mm=2^3;
%%%%%%%%%%%%%%
% End of the computation of the reference solution
%%%%%%%%%%%%%%
% Compute now 5 different approximations of the solution on
% different grids. n=1 coarser grid n=5 finest grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:5
    mm=mm*2;
    h=T/(mm);
    m=mm-1;
    hh=[hh h];
    xx=[h:h:T-h]';
    if (changeinval == 1 )
       theta0=0.7*ones(m,1)+sin(xx/2);
    else if (changeinval == 2)
        theta0=0.7*cos(xx);
        else
            theta0=0.7*cos(xx)+0.5*sin(xx);
        end
    end   
%
    theta=theta0;
    maxNit=100;
    figure(1)
    plot([0;xx;T],[alpha;theta0;beta], '--');
    grid
    hold on
    plot([0;xxx;T],[alpha;thetaex;beta]);
    drawnow
    title('solution and initial guess (dashed line)')
    hold off
%
% Newton iteration to solve G(theta)=0
%
    for k=1:maxNit
       thetaold=theta;
       J=JacG(theta);
       g=G(theta);
       theta=theta-J\g;
       % difference between two successive iterations as stopping criterion
       nr=norm(theta-thetaold)/norm(theta0)
       if (nr <tol )
           display('Number of Newton iterations to satisfy the convergence criterion')
           k
           break
       end    
    end
    % Number of Newton iterations to get convergence
    if (k==maxNit)
      display('Reached maximum number of iterations')
      maxNit
    end
%
% The grid used to compute the reference solution thetaex
% is finer than the grid used in each of the experiments of this for-loop
% we therefore map thetaex on the coarser grid to compare theta and thetaex
% and obtain an estimate of the error here called err()
%
  r=ceil((size(thetaex,1)+1)/mm);
  err(n)=norm((thetaex(r:r:end+1-r,1)-theta))*sqrt(h);
end
figure(3)
loglog(hh,err,'d-')
hold on
loglog(hh,hh.^2,'k')
err./(hh.^2)
xlabel('spatial step-size')
ylabel('2-norm of the error')
title('Log-log plot of the 2-norm of the error for a test BVP problem: order 2 discretization')
hold off