```

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;

```