%% 
% This code solves the nonlinear pendulum BVP 
% theta''+sin(theta) = 0 on an interval $t\in [0,T], 
% with boundary conditions theta(0) = alpha, theta(T) = beta. 

clear all, close all
clc

T=2*pi;
alpha=0.7; beta =0.7;

m=50; % Number of (interior) discretization points (equidistant)

h=T/(m+1);

t=linspace(0,T,m+2);
t=t';

%Theta = 0.7*cos(t)+0.5*sin(t);
%Theta=0.7*ones(m+2,1);
Theta=0.7+sin(t/2);

Theta=[alpha; Theta(2:end-1); beta];
error=1;
count=0;

while (error > 10^-1)
    
    plot(Theta,'--')
    G=[Theta(1:end-2)-2*Theta(2:end-1)+Theta(3:end)]+h^2*sin(Theta(2:end-1));

    B=zeros(m,3);
    B(:,1)=ones(m,1);
    B(:,2)=-2+h^2*cos(Theta(2:end-1));
    B(:,3)=ones(m,1);
    
    DG = spdiags(B,[-1 0 1],m,m);
    
    Delta=-DG\G;
    Thetanew=Theta+[0;Delta;0];
  
    error=norm(Theta-Thetanew);
    Theta=Thetanew;
    hold on
    count=count+1;
end

 plot(Theta,'-*')
count
