function twopoint1
%
% test program for modified exercise 7.3.1. of Van Loan, Introduction
% to Scientific Computing.  The modified method uses a compact
% discretization that is fourth-order accurate if p is constant.
% We print the solution at three points (x=0.25, x=0.50 and x=0.75) and
% also print the change in the solution at x=0.50 between successive
% computations at different mesh spacing.  From this one can judge the
% order of accuracy.  Should be asymptotically 2nd order for variable
% p and asymptotically 4th order for constant p.
%
  % Test run with the original function p of Ex. 7.3.1.
  fprintf(' original function p of 7.3.1\n')
  fprintf(' %6s %15s %15s %15s %9s\n', 'n',...
    'u(0.25)','u(0.50)','u(0.75)','du(0.50)');
  prev = 0;
  for k = 2:10
    n=1+2^k;
    x = linspace(0,1,n);
    u = CompactTwoPtBVP(n,0,1,'p','q','r');
    fprintf(' %6d %15.7e %15.7e %15.7e %10.2e\n', n,...
      u(1+(n-1)/4),u(1+(n-1)/2),u(1+3*(n-1)/4),u(1+(n-1)/2)-prev);
    prev = u(1+(n-1)/2);
  end
  fprintf('\n');
  % Test run with a modified function pmod
  fprintf(' modified function pmod, constant\n')
  fprintf(' %6s %15s %15s %15s %9s\n', 'n',...
    'u(0.25)','u(0.50)','u(0.75)','du(0.50)');
  prev = 0;
  for k = 2:10
    n=1+2^k;
    x = linspace(0,1,n);
    u = CompactTwoPtBVP(n,0,1,'pmod','q','r');
    fprintf(' %6d %15.7e %15.7e %15.7e %10.2e\n', n,...
      u(1+(n-1)/4),u(1+(n-1)/2),u(1+3*(n-1)/4),u(1+(n-1)/2)-prev);
    prev = u(1+(n-1)/2);
  end

  function [l,u] = TriDiLU(d,e,f)
%
% Pre:
%   d,e,f   n-vectors that define a tridiagonal matrix 
%           A = diag(e(2:n),-1) + diag(d) + diag(f(1:n-1),1).
%           A has an LU factorization.
% Post: 
%    l      n-vector with the property that 
%           L = eye + diag(l(2:n),-1)
%    u      n-vector with the property that 
%           U = diag(u) + diag(f(1:n-1),1).
%
   n = length(d); 
   l = zeros(n,1); 
   u = zeros(n,1);
   u(1) = d(1);
   for i=2:n
      l(i) = e(i)/u(i-1);
      u(i) = d(i) - l(i)*f(i-1);
   end

  function x = LBiDiSol(l,b)
%
% Pre:
%    l  n-vector that defines the lower bidiagonal matrix
%       L = I + diag(l(2:n),-1).
%    b  n-vector
%
% Post:
%    x  n-vector that solves Lx = b
%
   n = length(b); 
   x = zeros(n,1);
   x(1) = b(1);
   for i=2:n
      x(i) = b(i) - l(i)*x(i-1);
   end

  function x = UBiDiSol(u,f,b)
%
% Pre:
%  u,f  n-vectord that define the upper bidiagonal matrix
%       U = diag(u) + diag(f(1:n-1),1)
%    b  n-vector
%
% Post:
%    x  n-vector that solves Ux = b
%
   n = length(b); 
   x = zeros(n,1);
   x(n) = b(n)/u(n);
   for i=n-1:-1:1
      x(i) = (b(i) - f(i)*x(i+1))/u(i);
   end

  function uvals = CompactTwoPtBVP(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 midpoints; q and r at mesh points.
  p = zeros(1,n);
  q = zeros(1,n);
  r = zeros(1,n);
  for i = 1:n
    xi = a+(i-1)*h;
    if (i<n)
      p(i) = feval(pname,xi+h/2);
    end
    q(i) = feval(qname,xi);
    r(i) = feval(rname,xi);
  end
  % set up the matrix and right-hand-side for tridiagonal L*U solver.
  % d, f and e define the matrix, b defines the right-hand-side, all
  % in the notation of Van Loan, Section 6.2.1.
  d = zeros(1,n);
  e = zeros(1,n);
  f = zeros(1,n);
  b = zeros(1,n);
  % trivial equation at left boundary
  d(1) = 1;
  e(1) = 0;
  f(1) = 0;
  b(1) = 0;
  % interior equations
  for i = 2:n-1
    ??? d(i), e(i), f(i), b(i) ???
  end
  % trivial equation at right boundary
  d(n) = 1;
  e(n) = 0;
  f(n) = 0;
  b(n) = 0;
% factor and solve
  [tridl,tridu] = TriDiLU(d,e,f);
  y = LBiDiSol(tridl,b);
  uvals = UBiDiSol(tridu,f,y);

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

  function z = pmod(x)
%
% modified function p of exercise 7.3.1
%
  z = 2;

  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;