% This is a matlab m-file that
% solves the BVP below using finite differences
%
%   u'' + u = 0
%   u(0)=0
%   u(1)=beta
%
% (DMA, 5-28-07)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
a=0;
b=1;
beta=1;
ua=0;
ub=beta;
N=8;
dx=(b-a)/N;
dx2=dx*dx;
x=linspace(a,b,N+1)';    % creates N+1 discretized values of col. vector x
%
% define the matrix (no attempt to be efficient here ...)
%
A=zeros(N-1,N-1);
for i=1:N-1
    A(i,i)=-2+dx2;      % MAIN DIAGONAL
end
for i=1:N-2;
    A(i,i+1)=1;         % FIRST UPPER DIAGONAL
end
for i=2:N-1;
    A(i,i-1)=1;         % FIRST LOWER DIAGONAL
end
%
% set up the right hand side vector
%
b=zeros(N-1,1);
b(N-1,1)=-beta;
%
% solve A*u_int = b using Matlab's built in \ command
%
u_int=A\b;
%
% define the solution vector
%
u_sol=[ua;u_int;ub];
%
% define the known exact solution for comparison
%
uexact=sin(x)/sin(1);
%
% plot the approximate solution
%
plot(x,u_sol,'o-',x,uexact,'g--')
xlabel('x-axis');
ylabel('u-axis');
title('Finite Difference (o), Exact');

%
% define an error
%
norm(u_sol-uexact)/norm(uexact)

%
% error vs. N (shows second order accuracy)
%
%   N       dx            error norm
%   2       0.5          1.461101635926407e-003
%   4       0.25         4.145987091483337e-004
%   8       0.125        1.106235398883167e-004
%  16       0.0625       2.866937038673480e-005
%  32       0.03125      7.305453638474669e-006
%  64       0.015625     1.844418834717777e-006
% 128       0.0078125    4.634134925907454e-007