function solnMat = implicitFDHeatSolver(dt, dx, L, tmax, lBC, rBC, IC, gamma=1)
%function solnMat = implicitFDHeatSolver(dt, dx, L, tmax, lBC, rBC, IC, gamma=1)
%
%  Function uses a simple implicit (unconditionally 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)(:);
    rhs = uprev;
	rhs(1) += delta * lBC(j*dt);
	rhs(end) += delta * rBC(j*dt);
    ucurr = A \ rhs;
    solnMat = [solnMat; [lBC(j*dt) ucurr' rBC(j*dt)]];
  end

end