% [u,u_true,x,t] = beam_warming(t_0,t_f,M,N)
% 
% solves the advection equation u_t + c u_x = 0 on [0,2 pi]
%
% For initial data cos(x), we assume periodic initial data.  For
% step-function initial, we assume u_x = 0 at the ends.

function [u,u_true,x,t] = lax_friedrichs(t_0,t_f,M,N)

c = 1/2;

% define the mesh in space
dx = 2*pi/(N-1);
x = 0:dx:2*pi;
x = x';

% define the mesh in time
dt = (t_f-t_0)/M;
t = t_0:dt:t_f;

for i=1:N
  u(i,1) = f(x(i));
end 
u_true(:,1)=u(:,1);

for j=2:M+1
  u_old = u(:,j-1);

  for i = 3:N;
% first do the interior points: 
    u_new(i) = u_old(i)-c*dt/(2*dx)*(3*u_old(i)-4*u_old(i-1)+u_old(i-2)); 
    u_new(i) = u_new(i)+dt^2*c^2/(2*dx^2)*(u_old(i)-2*u_old(i-1)+u_old(i-2));
  end
% now do the endpoints u_new(1) and u_new(N).
% for periodic boundary conditions:
%  u_new(2) = u_old(2)-c*dt/(2*dx)*(3*u_old(2)-4*u_old(1)+u_old(N-1)); 
%  u_new(2) = u_new(2)+dt^2*c^2/(2*dx^2)*(u_old(1)-2*u_old(1)+u_old(N-1));
%  u_new(1) = u_old(1)-c*dt/(2*dx)*(3*u_old(1)-4*u_old(N-1)+u_old(N-2)); 
%  u_new(1) = u_new(1)+dt^2*c^2/(2*dx^2)*(u_old(1)-2*u_old(N-1)+u_old(N-2));
% for the step-function, we have u_old(2) = u_old(0) = u_old(-1), and
% u_old(N-1) = u_old(N+1).  If we kept track of ghost points at the
% 0 and N+1 indices.  Plugging this into the rule, we find:
  u_new(2)=u_old(1);
  u_new(1)=u_old(1);  

% now update the solution:
  u(:,j) = u_new';

% define the true solution: 
  t_diff = t(j)-t_0;
  for k=1:N
    u_true(k,j) = f(x(k)-c*t_diff);
  end
end