function [U Uexact]=implicit(dt,tf,dx,x);
% IMPLICIT Solves heat equation problem
%    u_t = u_xx, u(0,t)=u(1,t)=0, u(x,0)=x (1-x)
% by implicit method, exploiting sparseness.  Also returns exact solution.
% Example I:
%   >> dx=0.1; x=0:dx:1; dt=0.02; tf=0.5; [U, Ue]=implicit(dt,tf,dx,x);
%   >> plot(x,U,'o',x,Ue,'.'), xlabel x, legend('approx','exact')
% Example II:
%   >> dt=[1e-6 1e-5 5e-5 1e-4 5e-4 1e-3 5e-3 0.01 0.1];
%   >> tf=0.2;  dx=0.05;  x=0:dx:1;
%   >> for k=1:9, [U,Ue]=implicit(dt(k),tf,dx,x); err(k)=max(abs(U-Ue)); end
%   >> loglog(dt,err,'s'), xlabel('\Delta t')
%   >> ylabel(['error at t_f = ' num2str(tf)])
% ELB 2/24/05

J=length(x)-1;   nu=dt/(dx^2);
M=ceil(tf/dt); truetf=M*dt;
d=repmat(1+2*nu,1,J-1); % diagonal elements
s=repmat(-nu,1,J-2); % super- and sub-diagonal elements
A=sparse([1:J-1, 1:J-2, 2:J-1],[1:J-1, 2:J-1, 1:J-2],[d s s],J-1,J-1);

U=x.*(1-x);  % initial condition; U is J+1 row vector
for n=1:M
   b=U(2:J)'; % note  b   and  A\b  are column vectors;  U(2:J) is a row
   U(2:J)=(A\b)';
end
Uexact=uexact(x,truetf);

function y=uexact(x,t);
% Fourier series sum w/o loop!; x is row vector; t must be scalar
N=181;  % N such that sum is within 10^-6 of exact
mm=(1:2:2*N-1)';  pimm=pi*mm;
c=exp(-t*pimm.^2)./(mm.^3);
y=(8/pi^3)*sum(repmat(c,1,length(x)).*sin(pimm*x));