% This Matlab script solves the one-dimensional convection
% equation using a finite difference algorithm.  The 
% discretization uses central differences in space and forward
% Euler in time.  An inflow bc is set and at the outflow a
% one-sided (backwards) difference is used.  The initial condition
% is a Gaussian distribution.
%

clear all;

% Set-up grid
xL = -4;
xR =  4;
Nx = 401; % number of control volumes
x = linspace(xL,xR,Nx);

% Calculate cell size in control volumes (assumed equal)
dx = x(2) - x(1);

% Set velocity 
u = 1;

% Set final time
tfinal = 100;

% Set timestep
CFL = 0.1;
dt = CFL*dx/abs(u);

% Set initial condition
U = exp(-x.^2);
t = 0;

% Set bc state at left (assumes u>0)
UL = exp(-xL^2);

% Loop until t > tfinal
while (t < tfinal),

  % Copy old state vector
  U0 = U;
  
  % Inflow boundary
  U(1) = UL;
  
  % Interior nodes
  for i = 2:Nx-1,
    U(i) = U0(i) - dt*(U0(i+1)-U0(i-1))/(2*dx);
  end

  % Outflow boundary uses backward difference
  i = Nx;
  U(i) = U0(i) - dt*(U0(i)-U0(i-1))/(dx);
  
  % Increment time
  t = t + dt;

  % Plot current solution
  plot(x,U,'*');
  xlabel('x'); ylabel('U');
  title(sprintf('time = %f\n',t));
  axis([xL, xR, -0.5, 1.5]);
  grid on;
  drawnow;          
  
end