clear all
% A code that solves the linear wave equation 
% u_t + au_x = 0 by the upwind scheme
% using 
% 1. a cosine wave as initial data
% 2. step function

N = 101;
deltax = 1/(N-1);
x=[0:deltax:1]';
a=1;
% initial data:

% 1. cosine wave    
% u = 1 + cos(8*pi*x);


% if you want a cosine  function as initial data then cosinefun=1 and
% step=0!!!!
cosinefun=0;
if cosinefun
  u = 1 + cos(8*pi*x);
end

% if you want a step function as initial data then step=1 and cosinefun=0!!!!
step=1;
if step
  u=zeros(size(x));   
  I=find(x>0 & x<=0.5);
  u(I) = 1;
end 


% Courant number =a dt /  dx. Here you can change the courant number to see 
% how it effects the solution

CN = 1.0; 

dt=CN *deltax/a;
% nstp = 16;
t=0;
nstp=20; % 20 timesteps

for k = 1:nstp,
  t=t+dt; % Compute the time
  u = u-CN*(u-[u(N-1);u(1:end-1)]); % Update the solution by the upwind scheme
end;

% Exact solution in case of a step function
if step==1
  u_ex=zeros(size(x));
  I=find(x-a*t>0 & x-a*t<=0.5);
  u_ex(I) = 1;
elseif cosinefun==1
  u_ex = 1 + cos(8*pi*(x-a*t)); 
end

plot(x,u,'-+',x,u_ex,'--r')
legend('numerical solution','exact solution')

