  function u = BVP_FiniteDiff(a,ua,b,ub,q,r,n)
% a<b and ua and ub are scalars.
% q and r are handles to functions q(x) and r(x) defined on [a,b].
% n is a positive integer >=2. 

% Assume a<b and that n > 2. Let h = (b-a)/(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)
%   ------------    -  -----------            
%         h                h
% - ----------------------------------------  + q(i)u(i)  = r(i)
%                       h
% where 
%                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) + u(i)(2+h^2*q(i)) - u(i+1)  = h^2r(i)
%
%  We need to evaluate q and r...

x = linspace(a,b,n);
h = (b-a)/(n-1);
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) = 2+h^2*qVec(i+1);
    if i>1
        T(i,i-1) = -1;
    end
    if i<n-2
        T(i,i+1) = -1;
    end
    if i==1
       rhs(i) = h^2*rVec(i+1) + ua;
    elseif i==n-2
       rhs(i) = h^2*rVec(i+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;