%
%  lect28.m -- diffusion
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  parameters

D = 1;

N =   40;  L =   pi;  dx = 2*L/N;  dt = 1.0*(dx^2)/(2*D);
M =  400;  T = M*dt

%  set grids & coefficients

xg = 0:dx:2*L-dx;  x = [xg 2*L];
tg = 0:dt:T;
lam = D*dt/(dx^2);

%  IV initialization & plot

IVs  = exp(-(xg-pi).^2);
u = IVs;
usave = [u];

figure(1);  clf 
plot(xg,u,'r');  hold on;  plot(xg,u,'k.');
axis([0 2*pi -0.1 1.1]);

%  PDE solve with operation count

disp(' ')
disp(' simple diffusion scheme')
disp(' ')

for k=1:M 

% iteration
   u = u + lam*([u(2:end) u(1)] - 2*u + [u(end) u(1:end-1)]);

% interval plotting
   if (mod(k,M/4)==0);  
      usave = [usave ; u];
      plot(xg,u,'k');
   end

end

plot(xg,u,'b');  hold on;  plot(xg,u,'k.');

xlabel('\bf x-axis');  ylabel('\bf u-axis')
title('\bf diffusion')

figure(2)
waterfall(xg,[1:1:5],usave,ones(size(usave)))
axis([0 2*pi 1 5  0 1.1])
%axis equal

mean_u = mean(u)
max_minus_min = max(u)-min(u)