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;