function twopoint
% test program for exercise 7.3.1. of Charles Van Loan, Introduction to
% Scientific Computing
format long e;
hold on;
plotstyle = ['y-','m:','c-.','r--','g-','b:','k-.','y--','m-','c:',...
'r-.','g--','b-','k:','y-.','m--','c-','r:','g-.','b--','k-']
for k = 2:12
n=1+2^k
x = linspace(0,1,n);
u = TwoPtBVP(n,0,1,'p','q','r');
[u(1+(n-1)/4),u(1+(n-1)/2),u(1+3*(n-1)/4)]
plot(x,u,plotstyle(k))
pause(3)
end
hold off;

function [g,h] = CholTrid(d,e)
%
% Pre:
%    d,e    n-vectors with the property that the tridiagonal matrix
%           A = diag(d) + diag(e(2:n),-1) + diag(e(2:n),1) is positive
%           definite.
%
% Post:
%    g,h    n-vectors with the property that the lower bidiagonal
%           matrix G = diag(g) + diag(h(2:n),-1) satisfies A = GG^T.
%
n = length(d);
g = zeros(n,1);
h = zeros(n,1);
g(1) = sqrt(d(1));
for i=2:n
h(i) = e(i)/g(i-1);
g(i) = sqrt(d(i) - h(i)^2);
end;

function x = CholTridSol(g,h,b)
%
% Pre:
%    g    n-vector with nonzero entries
%    h    n-vector
%    b    n-vector
%
% Post:
%    x    column n-vector that solves GG'x = b where
%         G = diag(g) + diag(h(2:n),-1).
%
n = length(g);
y = zeros(n,1);
% Solve Gy = b
y(1) = b(1)/g(1);
for k=2:n
y(k) = (b(k) - h(k)*y(k-1))/g(k);
end
% Solve G'x = y
x = zeros(n,1);
x(n) = y(n)/g(n);
for k=n-1:-1:1
x(k) = (y(k) - h(k+1)*x(k+1))/g(k);
end;

function uvals = TwoPtBVP(n,a,b,pname,qname,rname)
%
% Pre:
%    a,b     reals with a<b.
%    n       integer >= 3.
%    pname   string that names a positive function defined on [a,b].
%    qname   string that names a positive function defined on [a,b].
%    rname   string that names a function defined on [a,b].
%
% Post:
%    uvals   column n-vector with the property that uvals(i) approximates
%            the solution to the two-point boundary value problem
%
%               -D[p(x)Du(x)] + q(x)u(x) = r(x)      a<=x<=b
%               u(a) = 0, u(b) = 0
%
%            at x = a+(i-1)*h where h = (b-a)/(n-1); 1<=i<=n
%
% compute mesh spacing
h = (b-a)/(n-1);
% evaluate p at mid-points, q and r at mesh points.
p = zeros(1,n-1);
q = zeros(1,n);
r = zeros(1,n);
for i = 1:n-1
xi = a+(i-1)*h;
p(i) = feval(pname,xi+h/2);
q(i) = feval(qname,xi);
r(i) = feval(rname,xi);
end
xi = b;
q(n) = feval(qname,xi);
r(n) = feval(rname,xi);
% set up the matrix and right-hand-side for CholTrid and CholTridSol.
% d and e define the matrix, b defines the right-hand-side.
d = zeros(1,n);
e = zeros(1,n);
b = zeros(1,n);
% trivial equation at left boundary
d(1) = 1;
e(2) = 0;
b(1) = 0;
% interior equations
for i = 2:n-1
d(i) = p(i-1)+p(i)+h*h*q(i);
if (i<n-1)
e(i+1) = -p(i);
end
b(i) = h*h*r(i);
end
% trivial equation at right boundary
d(n) = 1;
e(n) = 0;
b(n) = 0;
% factor and solve
[g,h] = CholTrid(d,e);
uvals = CholTridSol(g,h,b);

function z = p(x)
%
% function p of exercise 7.3.1
%
z = 2+x*cos(20*pi*x);

function z = q(x)
%
% function q of exercise 7.3.1
%
z = 20*sin(1/(0.1+x*x)^2);

function z = r(x)
%
% function r of exercise 7.3.1
%
z = 1000*(x-0.5)^3;