function steinerberger_test ( )

%*****************************************************************************80
%
%% steinerberger_test() tests steinerberger().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 June 2023
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../steinerberger' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'steinerberger_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test steinerberger().\n' );
%
%  Plots.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot f(n,x) for several values of n\n' );
  fprintf ( 1, '\n' );
  for n = [ 1, 5, 10, 25, 50, 100 ]
    x = linspace ( 0.0, 1.0, 5001 );
    fx = steinerberger_function ( n, x );
    clf ( );
    plot ( x, fx, 'linewidth', 2 );
    s = sprintf ( 'Steinerberger f(%d,x)', n );
    title ( s );
    grid ( 'on' );
    filename = sprintf ( 'steinerberger_%02d.png', n );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    close ( );
  end
%
%  Integral estimates.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Integrate f(n,x) for several values of n\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   n           exact        estimate      error\n' );
  fprintf ( 1, '\n' );
  xmin = 0.0;
  xmax = 1.0;
  for n = [ 1, 5, 10, 25, 50, 100 ]
    exact = 2.0 * harmonic_number ( n ) / pi;
    q = integral ( @(x) steinerberger_function ( n, x ), xmin, xmax );
    err = abs ( exact - q );
    fprintf ( 1, ' %3d  %14.6g  %14.6g  %14.6g\n', n, exact, q, err );
  end
%
%  Minimization.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Minimize f(n,x) for several values of n\n' );
  fprintf ( 1, '\n' );
% options = optimset ( 'Display', 'iter' );
  options = optimset ( );

  xmin = 0.21;
  xmax = 0.75;
  for n = [ 1, 5, 10, 25, 50, 100 ]
    [ x, fx, exitflag, output ] = fminbnd ( @(x) steinerberger_function ( ...
      n, x ), xmin, xmax, options );
    its = output.iterations;
    fprintf ( 1, ' %3d  %4d  %14.6g  %14.6g\n', n, its, x, fx );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'steinerberger_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  rmpath ( '../steinerberger' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 November 2022
%
%  Author:
%
%    John Burkardt
%
  t = now ( );
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
