```

function diffus
%
% Solve a diffusion equation.
%
% The equation is of the general form
%
%    (d/dt)u(t,x) + (d/dx)(-p(x)*(d/dx)u(t,x)) + q(x)*u(t,x) = r(x)
%
% on an interval a<=x<=b and 0<=t.
%
% The specific application here has p(x)=1, q(x)=1, r(x)=sin(x), and
% interval (a,b)=(0,pi).  The initial condition is u(0,x)=0.  The
% analytical solution is u(t,x)=0.5*sin(x)*(1-exp(-2*t)).
%
% The program integrates the equation from t=0 to t=1 using different
% mesh sizes and step numbers (defined by nx and nt).  The solution at
% the final time at x=pi/2 and the error at that point are printed.
%
fprintf(' %6s %6s %15s %10s\n', 'nx','nt','u(0.50)','err(0.50)');
t_end = 1;
exact = 0.5*(1-exp(-2*t_end));
for kx = 1:5
nx = 1+2^kx;
for kt = 2:10
nt = 2^kt;
dt = t_end/nt;
u = zeros(1,nx); % initial condition: u=0
for istep = 1:nt
u = StepDiffus(nx,0,pi,u,dt,'p','q','r');
end
fprintf(' %6d %6d %15.7e %10.2e\n', nx,nt,...
u(1+(nx-1)/2),u(1+(nx-1)/2)-exact);
end
fprintf('\n');
end
fprintf('\n');

function uvals = StepDiffus(n,a,b,u0,dt,pname,qname,rname)
%
% Pre:
%    n       integer >= 3.
%    a,b     reals with a<b.
%    u0      column n-vector that represents the solution at the start
%            of the present timestep.
%    dt      positive real value that is the timestep.
%    pname   string that names a positive function defined on [a,b].
%    qname   string that names a non-negative 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 diffusion equation at the end of the
%            present timestep.  The trapezoidal rule is used for time
%            integration, hence uvals is obtained by solving for u in:
%
%               u(x) + (dt/2)*(-D[p(x)Du(x)]+q(x)u(x))
%                 = u0(x) - (dt/2)*(-D[p(x)Du0(x)]+q(x)u0(x)) + dt*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
%
% A compact spatial discretization is employed that is fourth-order
% accurate in the case that p(x) is constant.
%
% 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 tridiagonal matrix for the spatial derivatives part
% of the problem.  d0, f0 and e0 define this tridiagonal matrix,
% in the notation of Van Loan, Section 6.2.1.
d0 = zeros(1,n);
e0 = zeros(1,n);
f0 = zeros(1,n);
% trivial equation at left boundary
d0(1) = 1;
e0(1) = 0;
f0(1) = 0;
% interior equations
for i = 2:n-1
d0(i) = (p(i-1)+p(i))+h*h*q(i)*10/12;
e0(i) = -p(i-1)+h*h*q(i-1)/12;
f0(i) = -p(i)+h*h*q(i+1)/12;
end
% trivial equation at right boundary
d0(n) = 1;
e0(n) = 0;
f0(n) = 0;
% 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 for diffusion equation.
%
z = 1;

function z = q(x)
%
% function q for diffusion equation.
%
z = 1;

function z = r(x)
%
% function r for diffusion equation.
%
z = sin(x);

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

```