 function uvals = TwoPtBVP(n,a,b,p,q,r,s)
 % 
 % a and b are reals with a<b and n is an integer with >= 3.
 % p, q, r, and s name functions defined on [a,b], the
 % first two of which are positive on the interval.
 % 
 % uvals is a column n-vector with the property that uvals(i) approximates
 % the solution to the 2-point boundary value problem
 %
 %                -D[p(x)Du(x)] + s(x)u'(x) + q(x)u(x) = r(x)    a<=x<=b
 %                 u(a) = u(b) = 0
 %
 % at x = a+(i-1)h where h = (b-a)/(n-1)
 
% Get the required evaluations of the functions p, s, q, and r
h = (b-a)/(n-1);
pvals = zeros(n,1); svals = zeros(n,1); qvals = zeros(n,1); rvals = zeros(n,1);
for i=2:n-1
    x = a+(i-1)*h;
    pvals(i) =  p(x-h/2);
    svals(i) =  s(x);
    qvals(i) =  q(x);
    rvals(i) =  r(x);
end
pvals(n) = p(b-h/2);


% First term discretization at x = a+(i-1)h ...
%        (pvals(i+1)(uvals(i+1)-uvals(i))/h  - pvals(i)(uvals(i)-uvals(i-1))/h )/h

% Second term discretization at x = a+(i-1)h ...
%     svals(i)( uvals(i+1)-uvals(i-1) )/(2h)

% Third term at x = a+(i-1)h ...
%     qvals(i)uvals(i)

% The linear equation at x = a+(i-1)h ...
%      e(i-1)uvals(i-1) + d(i)uvals(i) + f(i)uvals(i+1) = h^2 rvals(i)
% where
d = zeros(n,1); e = zeros(n,1); f = zeros(n,1);
for i=2:n-1
    e(i-1) = -pvals(i)   - (h/2)*svals(i);
    d(i)   =  pvals(i)   + pvals(i+1)    + h^2 * qvals(i);
    f(i+1) = -pvals(i+1) + (h/2)*svals(i);
end

% Here is what the system looks like for n = 7:

%     d(2)  f(3)   0    0    0    uvals(2)        rvals(2)
%     e(2)  d(3)  f(4)  0    0    uvals(3)        rvals(3)
%      0    e(3)  d(4) f(5)  0    uvals(4)  = h^2 rvals(4)
%      0     0    e(4) d(5) f(6)  uvals(5)        rvals(5)
%      0     0     0   e(5) d(6)  uvals(6)        rvals(6)

uvals = zeros(n,1);
uvals(2:n-1) = TriDiagSol(d(2:n-1),e(2:n-1),f(2:n-1),h^2*rvals(2:n-1));


function x = TriDiagSol(a,c,g,b)
% Solves Ax = b where A(i,i) = a(i), i=1:n, A(i+1,i) = c(i), i=1:n-1,
% and A(i,i+1) = g(i+1), i=1:n-1.

% Compute A = LU where L(i+1,i) = l(i), i=1:n-1, U(i,i) = u(i), i=1:n,
% and U(i,i+1) = g(i+1), i=1:n-1

n = length(b);
u(1) = a(1);
for i=1:n-1
    l(i) = c(i)/u(i);
    u(i+1) = a(i+1) - l(i)*g(i+1);
end

% Solve Ly = b
y(1) = b(1);
for i=2:n
    y(i) = b(i) - l(i-1)*y(i-1);
end


% Solve Ux = y
x(n) = y(n)/u(n);
for i=n-1:-1:1
    x(i) = (y(i) - g(i+1)*x(i+1))/u(i);
end




