nflux = { 'LaxFriedrichs', 'force', 'Richtmyer', ...
          'slic', 'GodunovCentered', 'ToroBillett', ...
          'Fromm', 'WarmingBeam', 'Waf' } ;
grid = [ 50 100 200 400 ];

for j=1:length(nflux)
  LinearAdvectionModel( 'setup',    -0.4 ) ;
  LinearAdvectionModel( 'set_test', 'smooth' ) ;
  LinearAdvectionModel( 'set_numerical_flux', nflux{j} ) ;
  LinearAdvectionModel( 'set_limiter', 'superbee' ) ;

  PRB . CFL       = 0.9*LinearAdvectionModel( 'get_cfl' ) ;
  PRB . FinalTime = LinearAdvectionModel( 'get_final_time' ) ;
  PRB . StepTime  = PRB . FinalTime / 100 ;

  SlopeLimiterModel( 'setup', 'superbee' ) ;
  % limiters:
  % 'minbee', 'superbee', 'vanleer', 'vanalbada', 'minmod', 'warmingbeam', 'laxwendroff', 'fromm'

  % call the solver
  for i=1:4
    PRB . Ncells = grid(i) ;
    [UFRAMES,TFRAMES] = Solver(PRB,@LinearAdvectionModel) ;
    %[UFRAMES,TFRAMES] = MUSCLSolver(PRB,@LinearAdvectionModel) ;
    [NR1(i),NR2(i),NRi(i)] = ComputeError(UFRAMES,TFRAMES,@LinearAdvectionModel);
  end
  disp( sprintf('\n\n%s', nflux{j}) ) ;
  for i=1:3
    disp( sprintf('norm1[%d,%d] = %f',grid(i+1),grid(i),log(NR1(i)/NR1(i+1))/log(2))) ;
    disp( sprintf('norm2[%d,%d] = %f',grid(i+1),grid(i),log(NR2(i)/NR2(i+1))/log(2))) ;
    disp( sprintf('normi[%d,%d] = %f',grid(i+1),grid(i),log(NRi(i)/NRi(i+1))/log(2))) ;
  end
end
