function exercise4 ( initial_condition )

%*****************************************************************************80
%
%% exercise4 compares upwind and Lax-Friedrichs solutions.
%
%  Discussion:
%
%    We consider BNI: the Burgers nonconservative inviscid equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function value = initial_condition ( x ): sets the initial condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4:\n' );
  fprintf ( 1, '  Compare upwind difference and Lax-Friedrichs solutions\n' );
  fprintf ( 1, '  for the Burgers equation.\n' );

  xmin = -1.0;
  xmax = +1.0;
  nx = 81;
  x = linspace ( xmin, xmax, nx );
  dx = ( xmax - xmin ) / ( nx - 1 );

  tmin = 0.0;
  tmax = 2.0;
  nt = 201;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );

  IM1 = [ nx-1, 1:nx-2 ];
  I   = [ 1:nx-1 ];
  IP1 = [ 2:nx-1, 1 ];

  for k = 1 : nt

    if ( k == 1 )
      un = initial_condition ( x );
      vn = initial_condition ( x );
%
%  Use both upwind (uo, un) and Lax-Friedrichs (vo, vn) schemes.
%
    else
      uo = un;
      un(I) = uo(I) - dt / dx * uo(I) .* ( uo(I) - uo(IM1) );
      un(nx) = un(1);

      vo = vn;
      vn(I) = ( vo(IP1) + vo(IM1) ) / 2.0 ...
        - ( dt / dx ) * vo(I) .* ( vo(IP1) - vo(IM1) ) / 2.0; 
      vn(nx) = vn(1);

    end
%
%  Plot the computed solutions.
%
    clf ( );
    plot ( x, un, 'b-', x, vn, 'r-', 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<--- x --->' );
    ylabel ( '<--- u(x,t) -->' );
    label = sprintf ( 'BNI Upwind vs Lax-Friedrichs: k = %3d, t = %4.1f', k, t(k) );
    title ( label );
    legend ( 'Upwind', 'Lax-Friedrichs' );
    axis ( [ xmin, xmax, 0.0, 1.1 ] );

    pause ( 0.1 );

  end

  filename = 'exercise4.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Last plot saved as "%s"\n', filename );
%
%  Report integral of difference.
%
  e = sum ( abs ( un(I) - vn(I) ) ) * dx;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For nx = %d, integral of solution difference = %g\n', nx, e );

  return
end

