function [x,t,u,Frames] = wave1d
% Transport and decay of pollutant

% Setup
L = 1.0;
T = 75;
d = 0.01;
%d = 0.0;
m = 150;
dt = T/m;
n =20;
dx = L/n;

v = 0.02;
u = zeros(n+1,m+1);

% Create mesh
x = (0:n)*dx;
t = (0:m)*dt;

% Set i.c. and b.c.s
u(:,1) = u0(x)';

u(1,:) = bcLeft;

alpha = v*dt/dx;

for j=1:m
    for i=2:n+1
        u(i,j+1) = alpha*u(i-1,j) + (1 - alpha -d*dt)*u(i,j);
    end
end

for l=1:m+1
    plot(x,u(:,l));
    axis([0 1 0 10]);
    Frames(l) = getframe;
end

% Some functions

function u0 = u0(x)
for i = 1:length(x)
   if x(i) < 0.5
       u0(i) = 5.0;
  else
       u0(i) = 0.0;
   end
end

function bcLeft = bcLeft
bcLeft = 0.0;
