clear all;
%% solution of heat equation with linear advection 
%
% Haberman, problem 10.4.3
%
%
% u_t = k*u_xx + c*u_x
% u(x,0) = f(x)
%
%  f(x) = 1, if abs(x) < a, 0 otherwise
%

%% parameters

c = 0.0; %   advecting velocity
k = 0.05; %   thermal diffusivity

a = 3; %   half-width of initial data

nX = 400;   % number of grid points for x domain
nT = 200;   % number of grid points for time domain

xmin = -50;
xmax = 20;
tmin = 0;
tmax = 200;

x1 = xmin:(xmax-xmin)/nX:xmax;
t1 = 0:tmax/nT:tmax;

[x,t] = meshgrid(x1,t1);

u = erf( ( x + a + c*t)./(2*sqrt(k*t)))/2 - erf( ( x-a + c*t)./(sqrt(4*k*t)))/2;

%%
figure(1);clf;
for j=1:5:tmax+1
   plot(x,u(j,:),'r-','LineWidth',2);
   set(gca,'FontSize',16);
   ylim([0,1.1]);
   titlestring = sprintf('u(x,t) at t = %d',(j-1)*tmax/nT);
   title(titlestring)
   xlabel('x');
   ylabel('u(x,t)');
   
   pause(0.01);
end

%%
figure(2);clf;
surf(x,t,u,'EdgeColor','none');
set(gca,'FontSize',16);
xlabel('x');
ylabel('t');
zlabel('u(x,t)')
titlestring = sprintf('advection-diffusion with c = %4.2f, k = %4.2f',c,k);
title(titlestring);

%%
figure(3);clf;
contourf(x,t,u,30,'LineColor','none');
set(gca,'FontSize',16);
xlabel('x');
ylabel('t');
h = colorbar;
set(h,'FontSize',16);
ylabel(h,'u(x,t)')

