%
% FTCS advection scheme
%
clear all

nx      =   201;		% number of grid points
D       =   40;			% width of domain
W0      =   -4;			% (backward) velocity
sigma   =   1;			% Gaussian profile width
Ampl    =   2;			% Gaussian profile amplitude
nt      =   200;		% number of timesteps
dt      =   1e-2;		% timestep
dx      =   D/(nx-1);		% grid step size
x       =   0:dx:D;		% 1D spatial domain

% Initial Gaussian T-profile
xc      =   20;			% center of Gaussian profile
T       =   Ampl*exp(-(x-xc).^2/sigma^2);

% Velocity
Vx      =   ones(1,nx)*W0;	% advection velocity
fprintf('courant = %g\n',...
	abs(W0)*dt/dx);		% Courant number, should be <1
cfac = dt/(2*dx); 

% Central finite difference discretization
for itime=1:nt    
    % central finite differences
    for ix=2:nx-1
      Tnew(ix) = T(ix) - Vx(ix)*cfac *(T(ix+1)-T(ix-1));
    end
    % BCs
    Tnew(1)     =   T(1);
    Tnew(nx)    =   T(nx);
    % Update Solution and Time Increment
    T           =   Tnew;
    time        =   itime*dt;
    % Analytical solution for this case
    T_anal      =   Ampl*exp(-(x-xc-time*W0).^2/sigma^2);
    % graphical plot
    figure(1),clf, plot(x,T,x,T_anal), legend('Numerical','Analytical')
    xlabel('x')
    ylabel('temperature')
    set(gca,'xlim',[0.0 30.0]);
    set(gca,'ylim',[-1.0 3.0]);
    % fprintf('%g of %g\n',itime,nt);
    drawnow
    % pause
end

