function flame_convergence_euler ( )

%*****************************************************************************80
%
%% flame_convergence_euler is a flame convergence study.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Cleve's Corner: Stiff Differential Equations,
%    MATLAB News and Notes,
%    May 2003, pages 12-13.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'flame_convergence_euler:\n' );
  fprintf ( 1, '  Show convergence of Euler method for the flame ODE.\n' );

  a = 0.0;
  b = 250.0;
  delta = 0.01;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot Euler solution, exact solution, and error.\n' );
  fprintf ( 1, '  Pause for 5 seconds after each plot.\n' );

  for n = [ 20, 40, 80, 160, 320, 640, 1280 ]
    h = ( b - a ) / n;
    [ t, y1 ] = flame_euler ( n, delta );
    y2 = flame_exact ( t, delta );
    plot ( t, y1, 'r-', t, y2, 'b-', t, abs(y1-y2), 'g-' );
    label = sprintf ( 'Forward Euler Approximation to Flame, N = %d\n', n );
    title ( label );
    legend ( 'Euler', 'Exact', 'Error' );
    grid ( 'on' );
    pause ( 5 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, ' n  h  e\n' );
  for n = [ 20, 40, 80, 160, 320, 640, 1280 ]
    h = ( b - a ) / n;
    [ t, y1 ] = flame_euler ( n, delta );
    y2 = flame_exact ( t, delta );
    e = rms ( y1 - y2 );
    fprintf ( 1, '  %4d  %12.4f  %12.4f\n', n, h, e );
  end

  fprintf ( 1, ' n  hold/h  eold/e\n' );
  for n = [ 20, 40, 80, 160, 320, 640, 1280 ]
    h = ( b - a ) / n;
    [ t, y1 ] = flame_euler ( n, delta );
    y2 = flame_exact ( t, delta );
    e = rms ( y1 - y2 );
    if ( 20 < n )
      fprintf ( 1, '  %4d  %12.4f  %12.4f\n', n, hold/h, eold/e );
    end
    hold = h;
    eold = e;
  end

  fprintf ( 1, ' n  log10(h)  log10(e)\n' );
  for n = [ 20, 40, 80, 160, 320, 640, 1280 ]
    h = ( b - a ) / n;
    [ t, y1 ] = flame_euler ( n, delta );
    y2 = flame_exact ( t, delta );
    e = rms ( y1 - y2 );
    fprintf ( 1, '  %4d  %12.4f  %12.4f\n', n, log10(h), log10(e) );
  end

  hvec = [];
  evec = [];
  for n = [ 20, 40, 80, 160, 320, 640, 1280 ]
    h = 1.0 / n;
    hvec = [ hvec, h ];
    [ t, y1 ] = flame_euler ( n, delta );
    y2 = flame_exact ( t, delta );
    e = rms ( y1 - y2 );
    evec = [ evec, e ];
  end

  plot ( log10(hvec), log10(evec), 'r-o', 'linewidth', 3 );
  axis ( 'equal' )
  grid ( 'on' );
  xlabel ( '<-- log10(h) -->' );
  ylabel ( '<-- log10(e) -->' );
  title ( 'Convergence rate' );
  filename = 'flame_convergence_euler.png';
  print ( '-dpng', filename );

  return
end

