% heateqn.m

clear all; close all;

% set various values
ell = 10 ;       % length of string
a=0.28  ;        % the left edge of the spike is at a*ell
d=0.3 ;          % the peak of the spike is at d*ell
b=0.32 ;         % the right edge of the spike is at b*ell
al=0.2 ;         % heat equation alpha (thermal diffusivity)
dx = 0.05 ;      % spacing in x space for plotting
dt = 0.02 ;      % time step for animation
tfinal = 5 ;     % final time for animation
nModes = 400 ;   % number of Fourier modes
delay = 0.0 ;    % extra time (in seconds) between animation frames 

% compute the Fourier coefficients for the initial temperature

beta = zeros(1, nModes) ;
for k = 1:nModes
     beta(k) = - 2*sin(k*pi*a)/( (d-a)*k*k*pi*pi ) ...
                + 2*(b-a)*sin(k*pi*d)/( (d-a)*(b-d)*k*k*pi*pi ) ...
                - 2*sin(k*pi*b)/( (b-d)*k*k*pi*pi ) ;
end

% compute the temporal decay rates

for k = 1:nModes
     Om(k) = (al*k*pi/ell)^2 ;
end

% run the time evolution

X = 0:dx:ell ;
set(gcf,'DoubleBuffer', 'on') ;
for t = 0:dt:tfinal             % for each fixed t
     B = beta .* exp(-t*Om) ;   % construct coeffs of sin(k*x*pi/ell) for time t
     U = zeros(1, length(X)) ;  % now construct u(x,t)
     for k = 1:nModes 
          U = U + B(k)*sin( k*pi*X/ell ) ;
     end
     plot(X,U) ;                % plot u(x,t)
     axis([0 ell -1 1]) ;
     title(sprintf('t = %4.1f', t)) ; 
     pause(delay) ;
end