  function ShowTwoPtBVP()
% Illustrates the function TwoPtBVP by apply it to the problem
%         du/dx^2 + u = 2x*sin(x) - 2*cos(x),   u(0) = u(pi) = 0
a = 0;  ua = 0;
b = pi; ub = 0;
n = 100;
u = TwoPtBVP(a,ua,b,ub,n,@p,@q,@r);
x = linspace(a,b,n)';
uExact = x.*sin(x);
err = max(abs(u-uExact));
plot(x,u,'r')
title(sprintf('h = %10.2e    err = %10.2e',(b-a)/(n-1),err),'fontsize',15)
shg

  function u = TwoPtBVP(a,ua,b,ub,n,p,q,r)
% Assume a<b and that n > 2. Let h = (b-1)/(n-1) and define x(i) = a+(i-1)h
% for i=1:n.
% u is a column n-vector with the property that u(i) approximates the
% solution to the following BVP at x(i):
%
%   -d/dx(p(x)du/dx) + q(x)u(x) = r(x)     u(a) = ua, u(b) = ub;
%

% Assume 2<=i<=n-1 and discretize at x(i):
%
%    u(i+1)-u(i)            u(i)-u(i-1)
%   ------------ p(i)   -  ----------- p(i-1)  + q(i)u(i)  = r(i)
%         h                       h
% - ----------------------------------------
%                       h
% where 
%                p(i) = p(x(i)+h/2) i=1:n-1
%                q(i) = q(x(i)),    i=2:n
%                r(i) = r(x(i)),    i=2:n
% 
%  Rearranging we get the following for i=2:n-1...
%
%       -u(i-1)p(i-1) + u(i)*(p(i-1)+p(i)+h^2*q(i)) - u(i+1)p(i)  = h^2r(i)
%
%  We need to evaluate p, q, and r...

x = linspace(a,b,n);
h = (b-a)/(n-1);

pVec = p(x(1:n-1)+h/2);
qVec = q(x);
rVec = r(x);

% Noting that u(1) = ua and u(n) = ub we assemble these equations in to 
% a system Tu(2:n) = rhs where

T = zeros(n-2,n-2);
rhs = zeros(n-2,1);
for i=1:n-2
    T(i,i) = pVec(i)+pVec(i+1)+h^2*qVec(i+1);
    if i>1
        T(i,i-1) = -pVec(i);
    end
    if i<n-2
        T(i,i+1) = -pVec(i+1);
    end
    if i==1
       rhs(i) = h^2*rVec(i+1) + pVec(i)*ua;
    elseif i==n-2
       rhs(i) = h^2*rVec(i+1) + pVec(n-1)*ub;
    else
       rhs(i) = h^2*rVec(i+1);
    end
end
u = zeros(n,1);
u(1) = ua;
u(n) = ub;
u(2:n-1) = T\rhs;

function y = p(x)
y = ones(length(x),1);

function y = q(x)
y = ones(length(x),1);


function y = r(x)
y = 2*x.*sin(x) - 2*cos(x);


        
      
        



