function solnMat = fdHeatSolver(dt, dx, L, tmax, lBC, rBC, IC, gamma=1)
%function solnMat = fdHeatSolver(dt, dx, L, tmax, lBC, rBC, IC, gamma=1)
%
%  Function uses a simple explicit (conditionally stable) finite
%  difference algorithm to solve the one-dimensional homogeneous
%  constant-coefficient heat equation on a bounded interval with
%  Dirichlet boundary conditions.
%
%  INPUTS:
%    dt     The temporal step size---should divide evenly into tmax.
%    dx     The spatial step size---should divide evenly into L.
%    L      The right end of the spatial interval [0, L] for problem.
%    tmax   The final time for which a solution is desired.
%    lBC    Function handle which gives the left boundary condition at x = 0.
%           This, along with other functions supplied as arguments, needs
%           to be one which can evaluate a vector of t-values.
%    rBC    Function handle which gives the right boundary condition at x = L.
%    IC     Function handle which gives the initial condition at t = 0.
%    gamma  Diffusivity coefficient (constant); defaults to 1
%
%  OUTPUT:
%  solnMat  A matrix whose first row is the discrete solution at time
%           t = 0, 2nd row the solution at time t = dt, 3rd row the
%           solution at time t = 2*dt, etc.

delta = gamma*dt/dx^2;
N = L / dx;
A = diag((1-2*delta)*ones(N-1,1)) + diag(delta*ones(N-2,1), 1) ...
     + diag(delta*ones(N-2,1), -1);

solnMat = IC(0:dx:L)(:)';
for j=1:(tmax/dt)
  uprev = solnMat(j, 2:N)(:);
  ucurr = A * uprev;
  ucurr(1) += delta * lBC((j-1)*dt);
  ucurr(N - 1) += delta * rBC((j-1)*dt);
  solnMat = [solnMat; [lBC(j*dt) ucurr' rBC(j*dt)]];
end

end