%   _                    _____     _          _      _      _         
%  | |    __ ___  __    |  ___| __(_) ___  __| |_ __(_) ___| |__  ___ 
%  | |   / _` \ \/ /____| |_ | '__| |/ _ \/ _` | '__| |/ __| '_ \/ __|
%  | |__| (_| |>  <_____|  _|| |  | |  __/ (_| | |  | | (__| | | \__ \
%  |_____\__,_/_/\_\    |_|  |_|  |_|\___|\__,_|_|  |_|\___|_| |_|___/
%
% Compute the Lax-Friedrichs flux
%
function Flux = LaxFriedrichs( MODEL, DT, DX, US )
  FS = feval( MODEL, 'flux', US ) ;
  %
  UL = US(:,2:end-2) ;
  UR = US(:,3:end-1) ;
  FL = FS(:,2:end-2) ;
  FR = FS(:,3:end-1) ;
  %
  Flux = 0.5 .* ( FL + FR + (DT/DX) * (UL-UR) ) ;
end
