function [t,x,U] = heat1 (T,a,m,n,beta,c,f,g)
%----------------------------------------------------------------------
% Usage:       [t,x,U] = heat1 (T,a,m,n,beta,c,f,g)
%
% Description: Use the Crank-Nicolson method to solve the following
%              partial differential equation called the one-dimensional
%              heat equation (also the diffusion equation):
%
%                  (d/dt)u(t,x) = beta*(d/dx)*(d/dx)u(t,x)
%
%              The heat equation is to be solved over the region:
%
%                  0 <= t <= T,
%                  0 < x < a
%
%              subject to the following initial and boundary conditions:
%
%                  u(0,x) = f(x)
%                  (1 - c(1))*u(t,0) - c(1)*(d/dt)u(t,0) = g(t,1)
%                  (1 - c(2))*u(t,a) + c(2)*(d/dt)u(t,a) = g(t,2)
%  
%              where c(1) and c(2) are in the range [0,1].
%
% Inputs:      T    = final solution time (T > 0)
%              a    = upper boundary of x (a > 0)
%              m    = number steps in t (m >= 1)
%              n    = number of steps in x (n >= 1)
%              beta = heat equation coefficient (beta > 0)
%              c    = 1 by 2 vector specifying boundary conditions
%                     0 = Dirichlet, 1 = Neumann (see above).
%              f    = string containing name of user-supplied 
%                     initial condition function.
%              g    = string containing name of user-supplied 
%                     boundary condtion function 
%                       
%                     The functions f and g are of the form:
%
%                     function y = f(x)
%                     function y = g(t,k)
%
%                     When f is called, it must return the initial value
%                     u(0,x).  When g is called, it must return the
%                     boundary value g(t,k). The functions f and g are
%                     optional in the sense that either can be replaced
%                     by the empty string ''. When this is done, the 
%                     corresponding function is assumed to be zero
%                     which means that a user-supplied function is
%                     NOT required in this case.
%
% Outputs:     t = (m+1) by 1 vector of t grid points:
%                  t(k) = (k-1)*T/m for 1 <= k <= m+1.
%              x = n by 1 vector of x grid points:
%                  x(j) = ja/(n+1) for 1 <= j <= n. 
%              U = (m+1) by n matrix containing the estimated solution
%                  with U(k,j) = u(t(k),x(j)).
%
% Notes:       1.  The function heat1 uses the functions trifac and trisub
%                  to solve a tridiagonal linear algebraic system at each
%                   time step.
%
%              2.  The Crank-Nicolson method is a stable method that is
%                  second-order accurate in both time and space. 
%----------------------------------------------------------------------
   
% Initialize 

   T = args (T,eps,T,1,'heat1');
   a = args (a,eps,a,2,'heat1');
   m = args (m,1,m,3,'heat1');
   n = args (n,1,n,4,'heat1');
   beta = args (beta,eps,beta,5,'heat1');
   A = zeros (3,n);
   Q = zeros (3,n);
   b = zeros (n,1);
   dt = T/m;
   dx = a/(n+1);
   t = [0 : dt : T]';
   x = [dx : dx : n*dx]';
   U = zeros (m+1,n);
   v = zeros (n,1);
   gamma = beta*dt/(dx*dx);   
   funf = getfun (f);
   fung = getfun (g);

% Do row 1

   if funf
      for j = 1 : n
         U(1,j) = feval(f,j*dx);
      end
   end
         
% Compute coefficient matrix of u(k+1) 

   b1 = c(1) + (1 - c(1))*dx;
   b2 = c(2) + (1 - c(2))*dx;
   g0 = 2*(1 + gamma);
   for k = 1 : n
      A(1,k) = -gamma;
      A(2,k) =  g0;
      A(3,k) = -gamma;
   end
   A(2,1) = A(2,1) - gamma*c(1)/b1;
   A(2,n) = A(2,n) - gamma*c(2)/b2;

   [Q,delta] = trifac (A); 
   if abs(delta) < eps/100
      disp ('Tridiagonal coefficient matrix in heat1 is singular')
   end

% Solve the PDE 

   g0 = 2*(1 - gamma);
   for k = 1 : m
      b(1) = (g0 + (c(1)/b1)*gamma)*U(k,1) + gamma*U(k,2);
      for j = 2 : n-1
         b(j) = gamma*U(k,j-1) + g0*U(k,j) + gamma*U(k,j+1);
      end
      b(n) = gamma*U(k,n-1) + (g0 + (c(2)/b2)*gamma)*U(k,n);
      if fung
         b(1) = b(1) + (gamma*dx/b1)*(feval(g,k*dt,1) + ... 
                feval(g,(k+1)*dt,1));
         b(n) = b(n) + (gamma*dx/b2)*(feval(g,k*dt,2) + ...
                feval(g,(k+1)*dt,2));
      end
      [v,delta] = trisub(Q,b);
      U(k+1,1:n) = v';
   end
%----------------------------------------------------------------------