%
% twopbvp_1.m 
% second order finite difference method for the two point b.v.p.
%   -u''(x) + q(x) u(x) = f(x),   u(ax) = u(bx) = 0.
% (Uses 3-pt differences on a uniform grid.)
%
clear all;
close all;
%
% N = # of intervals,
% ax,bx = the endpoints of the interval [a,b]
%
N = 10;
ax = 0;
bx = 1;
h = (bx-ax)/N;
x = linspace(ax+h,bx-h,N-1);
%
% Define functions using function handles
%
ff = @(x) sin(2*pi*x) .* (cos(pi*x) + 1 + 4*pi*pi);  % right hand side function
qq = @(x) cos(pi*x) + 1;                             % function q(x)
uexact = @(x) sin(2*pi*x);                           % exact soln
%
% Construct the coefficient matrix
%
d = 2 + h * h * qq(x);
dl = -ones(1,N-2);
A = diag(d) + diag(dl,-1) + diag(dl,1);
%
% Compute the right hand side vector
%
b = (h * h * ff(x))';
%
minTime = Inf;
tstart = tic;
%
u = A\b;                          % Solve the linear system A*x = b
%
telapsed = toc(tstart);           % Estimate elapsed time for the 
minTime = min(telapsed,minTime);  % solving the linear system
%
fprintf(1, '\n');
fprintf(1, 'time elapsed = %9.5e\n', minTime);
%
% Plot the approximate vs. the exact solution
%
xx = linspace(ax,bx,101);
plot(x,u,'-o',xx,uexact(xx),'r');
legend('approx.','exact');
grid on;
title('Numerical soln of 2-point b.v.p. -(d^2u)/(dx^2) + qu = f, u(a)=u(b)=0')
xlabel('x')
%
ex = (uexact(x))';
err = norm(u-ex,inf)/norm(ex,inf);     % Compute the relative max error
%
fprintf(1, '\n');
fprintf(1, 'relative max error = %12.8f\n', err);
