function U = implicitheat(J,N,tf)
% IMPLICITHEAT Solves u_t = u_xx by simplest implicit (\theta = 1)
% method.
% Examples:
%    >> implicitheat(10,10,0.01)
%    >> implicitheat(1000,1000,0.1);

% ELB 2/26/07

dx = 1/J;  mu = (tf/N) / (dx*dx)
x = 0:dx:1;  
U = x .* (1-x);

e = ones(J-1,1);
A = spdiags([-mu*e, (1+2*mu)*e, -mu*e], -1:1, J-1, J-1);
%full(A)

for n=1:N
    %   system is   A * U(2:J) = Uold(2:J)
    U(2:J) = (A \ U(2:J)')';
    %same as:   U(2:J) = (U(2:J) / A)';
end

plot(x,x.*(1-x), x,U,'--');