function c_new = Lax_Wendroff(c_old,cfl,Lim)
% executes one limited-Lax-Wendroff time step

% compute fluxes  flx is assumed staggered dx/2 to the right of c_old
flx = cfl*c_old + cfl*(1-cfl)/2*(circshift(c_old,-1)- ...
        circshift(c_old,0) ).*Lim;

% update solution
c_new = c_old - flx + circshift(flx,1);



