function [uu,tau]= heat(fc,nstep,tau)
% Fourier solution of 1-D heat equation
%   input parameters
%          fc       row vector of fourier coefficients of f(x)
%          nstep    number of time steps
%          tau      final time
ndiv=256;  h=1/ndiv;  x=0:h:1; % set spatial resolution
dt=tau/nstep; tau=0; % set temporal resolution
% define and plot initial condition
ell=length(fc); j=1:ell; X=sin(j.*pi.*x'); ww=fc.*X;
uu=sum(ww'); stepsolplot(tau,x,uu), shg
% timestep loop
for step=1:nstep
tau=tau+dt;  T=exp(-j.*j.*pi.*pi.*tau); ww=fc.*X.*T;
uu=sum(ww'); stepsolplot(tau,x,uu)
end
                                             
function stepsolplot(tau,x,u)
if tau==0, % plot initial condition
plot(x,u,'-b'), axis('square'), title('initial profile'),
xlabel('x'), ylabel('temperature'), shg, pause(2)
else % plot solution
hold on, plot(x,u,'-b'),
title(['temperature: time is ',num2str(tau)])
pause(0.5), hold off
end
return
