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