function [x, t, u] = heateq(f,g1,g2,L,T,m,n,c)
%heateq Solve and plot the solution of the heat equation using
%the standard finite difference method with accuracy O(h^2) for the
%spacial step and O(k) for the time step.
%
%   [x t u] = heateq(f,g1,g2,L,T,m,c) approximates the solution
%   to u_t = c u_xx for 0<=x<=L, 0<=t<=T with
%   boundary conditions: u(0, t) = g1(t); u(L,t) = g2(t),
%   initial Condition:   u(x, 0) = f(x)
%   Input parameters:
%     f -- function for initial condition
%     g1,g2 -- functions for boundary conditions
%     L -- length of rod
%     T -- length of time interval
%     m -- number of subintervals for x
%     n -- number of subintervals for t
%     c -- rate constant in equation
%   Output:
%     t -- vector of time points
%     x -- vector of x points
%     u -- matrix of the solution, u(i,j) ~= u(x(i),t(j))

 % spacial increment
 h = L/m;
 % time increment
 k = T/n;
 % compute constant
 r = c*k/h^2;
 % discretize in space and time
 x = (0:h:L)';  
 t = (0:k:T)';
 % compute the three diagonals
 ld = -r*ones(m-2,1);
 d = (1+2*r)*ones(m-1,1);
 ud = -r*ones(m-2,1);

 % create the matrix for u
 u = zeros(m+1,n+1);
 % create the rhs for b
 b = zeros(m-1,1);

 % evaluate initial conditions
 u(:,1) = feval(f,x);

 % evaluate boundary conditions
 u(1,:) = feval(g1,t);
 u(m+1,:) = feval(g2,t);

% find the solution at the remaining time steps
%  for j = 1:n
% 	  b(1) = r*u(j+1,1);
% 	  b(m-1) = r*u(j+1,m+1);
% 	  u(j+1,2:m) = thomas(ld, d, ud, b);
%  end
for j = 1:n
   i = 2;
   b(1) = u(i,j) + r*u(i-1,j+1);
   for i = 3:m-1
      b(i-1) = u(i,j);
   end
   b(m-1) = u(m,n) + r*u(m+1,n+1);
	u(2:m,j+1) = thomas(ld, d, ud, b);
end

 mesh(x,t,u');
 strtitle = sprintf('Plot of the solution,\n0 <= x <= %.2f, 0 <= t <= %.2f',...
							L, T);
 title(strtitle);
 xlabel('x-axis');
 ylabel('t-axis');
 zlabel('Temperature');
