function stiff_convergence_backward ( )

%*****************************************************************************80
%
%% stiff_convergence_backward looks at backward euler convergence for stiff ODE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'stiff_convergence_backward:\n' );

  a = 0.0;
  b = 1.0;

  fprintf ( 1, ' n  h  e\n' );
  for n = [ 20, 40, 80, 160, 320, 640 ]
    h = ( b - a ) / n;
    [ t, y1 ] = stiff_euler_backward ( n );
    y2 = stiff_exact ( t );
    e = rms ( y1 - y2 );
    fprintf ( 1, '  %3d  %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 ]
    h = ( b - a ) / n;
    [ t, y1 ] = stiff_euler_backward ( n );
    y2 = stiff_exact ( t );
    e = rms ( y1 - y2 );
    if ( 20 < n )
      fprintf ( 1, '  %3d  %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 ]
    h = ( b - a ) / n;
    [ t, y1 ] = stiff_euler_backward ( n );
    y2 = stiff_exact ( t );
    e = rms ( y1 - y2 );
    fprintf ( 1, '  %3d  %12.4f  %12.4f\n', n, log10(h), log10(e) );
  end

  hvec = [];
  evec = [];
  for n = [ 20, 40, 80, 160, 320, 640 ]
    h = 1.0 / n;
    hvec = [ hvec, h ];
    [ t, y1 ] = stiff_euler_backward ( n );
    y2 = stiff_exact ( t );
    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 = 'stiff_convergence_backward.png';
  print ( '-dpng', filename );

 return
end

