function zoomin_all_test ( a, b, c, abserr, kmax, f, ipoly, mult, nder, nsub )

%*****************************************************************************80
%
%% zoomin_all_test() calls all the zero finders for a given problem.
%
%  Discussion:
%
%    The original code was written by Harold Deiss in BASIC.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A, B, C, three estimates for a root of the
%    function.  The bisection routines will only be called if the function
%    evaluated at A, B and C includes both positive and negative values.
%
%    real ABSERR, an error tolerance.
%
%    function value = F(x,ider), the name of the routine that
%    evaluates the function or its derivatives.
%
%    integer IPOLY, is 0 if the function is not known to be
%    a polynomial, or is equal to the polynomial degree otherwise.
%
%    integer KMAX, the maximum number of iterations allowed.
%
%    integer MULT, the multiplicity of the root being 
%    sought.  If the multiplicity of the root is not known, set MULT to 1.
%
%    integer NDER, the highest order derivative that can be
%    computed by the user function.  NDER = 0 means only the function
%    value itself is available, for instance.  NDER must be at least 0,
%    and no algorithm needs a derivative higher than 3.
%
%    integer NSUB, the number of substeps to take.  A few
%    algorithms include a "sub-iteration".  For instance, the modified
%    Newton method evaluates the derivative function only before
%    steps 1, NSUB+1, 2*NSUB+1 and so on, and uses that derivative
%    value for NSUB iterations in a row.  This option is useful when
%    the derivative evaluation is expensive.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'zoomin_all_test():\n' );
  fprintf ( 1, '  A compilation of scalar zero finders,\n' );
  fprintf ( 1, '  based on the work of Joseph Traub.\n' );

  fa = f ( a, 0 );
  fb = f ( b, 0 );
  fc = f ( c, 0 );
%
%  Rearrange the data, if necessary, so that the pair (A,FA) has
%  the smallest function value of the three sets of data.
%
  if ( abs ( fb ) < abs ( fa ) )
    t = a;
    a = b;
    b = t;
    t = fa;
    fa = fb;
    fb = t;
  end

  if ( abs ( fc ) < abs ( fa ) )
    t = a;
    a = c;
    c = t;
    t = fa;
    fa = fc;
    fc = t;
  end
%
%  If necessary, switch B and C, so that FB is of opposite sign to FA,
%  if possible.
%
  sa = sign ( fa );
  sb = sign ( fb );
  sc = sign ( fc );

  if ( sa ~= sb )
    change = true;
  elseif ( sa ~= sc )
    t = b;
    b = c;
    c = t;
    t = fb;
    fb = fc;
    fc = t;
    change = true;
  else
    change = false;
  end
%
%  Check other input.
%
  if ( mult < 1 )
    mult = 1;
  end

  if ( nder < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'zoomin_all_test(): Fatal error!\n' );
    fprintf ( 1, '  The input quantity NDER = %d\n', nder );
    fprintf ( 1, '  But NDER must be at least 0.\n' );
    return
  end

  if ( nsub < 1 )
    nsub = 1;
  end
%
%  Report start up conditions.
%
  fprintf ( 1, '  1 point formulas use:\n' );
  fprintf ( 1, '    x1 = %g\n', a );
  fprintf ( 1, '    fx1= %g\n', fa );
  fprintf ( 1, '  2 point formulas add:\n' );
  fprintf ( 1, '    x2 = %g\n', b );
  fprintf ( 1, '    fx2= %g\n', fb );
  fprintf ( 1, '  3 point formulas add:\n' );
  fprintf ( 1, '    x3 = %g\n', c );
  fprintf ( 1, '    fx3= %g\n', fc );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  User estimated multiplicity = %d\n', mult );
  if ( 0 <= ipoly )
    fprintf ( 1, '  Polynomial degree = %d\n', ipoly );
  else
    fprintf ( 1, '  The function is not known to be polynomial.\n' );
  end
  fprintf ( 1, '  Highest derivative supplied = %d\n', nder );
  fprintf ( 1, '  Error tolerance = %g\n', abserr );
  fprintf ( 1, '  Maximum number of steps = %d\n', kmax );
  fprintf ( 1, '  Newton method substep parameter  = %d\n', nsub );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Technique                           Root        Steps' );
  fprintf ( 1, '   Error     Multiplicity\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  1.  One point iteration functions with memory:\n' );
  fprintf ( 1, '\n' );
%
%  Secant algorithm.
%
  name = 'Secant';
  x = a;
  x1 = b;
  [ x, x1, ierror, k ] = secant ( x, x1, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Extended secant algorithm.
%
  name = 'Extended secant';
  x = a;
  x1 = b;
  x2 = c;
  [ x, x1, x2, ierror, k ] = secant_extended ( x, x1, x2, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Capital Phi(2,1).
%
  name = 'Capital Phi(2,1)';
  x = a;
  x1 = b;
  x2 = c;
  [ x, x1, x2, ierror, k ] = cap_phi_21 ( x, x1, x2, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Muller's algorithm.
%
  name = 'R8_Muller';
  x = a;
  x1 = b;
  x2 = c;
  [ x, x1, x2, ierror, k ] = r8_muller ( x, x1, x2, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Perp E(2,1).
%
  if ( 1 <= nder )

    name = 'Perp E(2,1)';
    x = a;
    x1 = b;
    x2 = c;
    [ x, x1, x2, ierror, k ] = perp_e_21 ( x, x1, x2, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Star E(2,1).
%
  name = 'Star E(2,1)';
  x = a;
  x1 = b;
  x2 = c;
  [ x, x1, x2, ierror, k ] = star_e21 ( x, x1, x2, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Finite difference Halley's method.
%
  name = 'Finite difference Halley';
  x = a;
  x1 = b;
  x2 = c;
  [ x, x1, x2, ierror, k ] = halley2 ( x, x1, x2, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Phi(1,2).
%
  if ( 1 <= nder )

    name = 'Phi(1,2)';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = phi_12 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Perp E(1,2).
%
  if ( 1 <= nder )

    name = 'Perp E(1,2)';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = perp_e_12 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Star E(1,2).
%
  if ( 1 <= nder )

    name = 'Star E(1,2)';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = star_e12 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Dagger E(1,2).
%
  if ( 1 <= nder )

    name = 'Dagger E(1,2)';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = dagger_e12 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  2. One point iteration functions.\n' );
  fprintf ( 1, '\n' );
%
%  Newton method.
%
  if ( 1 <= nder )

    name = 'Newton';
    x = a;
    [ x, ierror, k ] = newton ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Steffenson's method.
%
  name = 'Steffenson';
  x = a;
  [ x, ierror, k ] = steffenson ( x, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Stirling's method.
%
  if ( 1 <= nder )

    name = 'Stirling';
    x = a;
    [ x, ierror, k ] = stirling ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Midpoint method.
%
  if ( 1 <= nder )

    name = 'midpoint';
    x = a;
    [ x, ierror, k ] = midpoint ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Traub-Ostrowski method.
%
  name = 'Traub-Ostrowski';
  x = a;
  [ x, ierror, k ] = traub_ostrowski ( x, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Chebyshev method.
%
  if ( 2 <= nder )

    name = 'Chebyshev';
    x = a;
    [ x, ierror, k ] = chebyshev ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Halley Super method.
%
  if ( 2 <= nder )

    name = 'Halley Super';
    x = a;
    [ x, ierror, k ] = halley_super ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Whittaker convex acceleration method.
%
  if ( 2 <= nder )

    name = 'Whittaker';
    x = a;
    [ x, ierror, k ] = whittaker ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Whittaker double convex acceleration method.
%
  if ( 2 <= nder )

    name = 'Whittaker2';
    x = a;
    [ x, ierror, k ] = whittaker2 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  E3.
%
  if ( 2 <= nder )

    name = 'E3';
    x = a;
    [ x, ierror, k ] = e3 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  E4.
%
  if ( 3 <= nder )

    name = 'E4';
    x = a;
    [ x, ierror, k ] = e4 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Halley's method.
%
  if ( 2 <= nder )

    name = 'Halley';
    x = a;
    [ x, ierror, k ] = halley1 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  PSI(2,1).
%
  if ( 3 <= nder )

    name = 'Psi(2,1)';
    x = a;
    [ x, ierror, k ] = psi_21 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  PSI(1,2).
%
  if ( 3 <= nder )

    name = 'Psi(1,2)';
    x = a;
    [ x, ierror, k ] = psi_i2 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Capital PHI(0,3).
%
  if ( 2 <= nder )

    name = 'Capital Phi(0,3)';
    x = a;
    [ x, ierror, k ] = cap_phi_03 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Reduced capital PHI(0,4).
%
  if ( 3 <= nder )

    name = 'Reduced Capital Phi(0,4)';
    x = a;
    [ x, ierror, k ] = red_cap_phi_04 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Functions from Hansen and Patrick.
%
%  The Ostrowski square root formula is the Hansen family member with BETA = 0.
%
  if ( 2 <= nder )
    name = 'Ostrowski square root';
    x = a;
    [ x, ierror, k ] = ostrowski_sqrt ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  The Euler method.
%
  if ( 2 <= nder )

    name = 'Euler';
    x = a;
    [ x, ierror, k ] = euler ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  The Laguerre method.
%
  if ( 2 <= nder && 2 <= ipoly )

    name = 'Laguerre';
    x0 = a;
    [ x, ierror, k ] = laguerre ( x0, ipoly, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  3.  Multipoint iteration functions.\n' );
  fprintf ( 1, '\n' );
%
%  First function in Traub, page 236.
%
  if ( 1 <= nder )

    name = 'Traub first';
    x = a;
    [ x, ierror, k ] = traub1 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  A family of iteration functions including
%  forms 2, 12, and 13 of Traub, pages 236 - 237.
%
  if ( 1 <= nder )

    name = 'Traub second';
    e = 0.5;
    d = 1.0;
    x = a;
    [ x, ierror, k ] = t_family1 ( x, e, d, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'Traub twelfth';
    e = 0.25;
    d = 2.0 / 3.0;
    x = a;
    [ x, ierror, k ] = t_family1 ( x, e, d, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'Traub thirteenth';
    e = 5.0 / 12.0;
    d = 6.0 / 7.0;
    x = a;
    [ x, ierror, k ] = t_family1 ( x, e, d, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Modified Newton method.
%
  if ( 1 <= nder && 1 <= nsub )

    name = sprintf ( 'Modified Newton, NSUB = %d', nsub );
    x = a;
    [ x, ierror, k ] = newton_mod ( x, nsub, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Fourth function in Traub.
%
  if ( 1 <= nder && 1 <= nsub )

    name = sprintf ( 'Traub fourth, NSUB = %d', nsub );
    x = a;
    [ x, ierror, k ] = traub4 ( x, nsub, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Newton - secant.
%
  if ( 1 <= nder )

    name = 'Newton - secant';
    x = a;
    [ x, ierror, k ] = newton_sec ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  A family of iteration functions including
%  forms 6 and 7 of Traub, pages 236 - 237.
%
  if ( 1 <= nder )

    name = 'Traub sixth';
    e = 2.0;
    ff = 3.0;
    g = 1.0;
    h = 1.0;
    x = a;
    [ x, ierror, k ] = t_family2 ( x, e, ff, g, h, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'Traub seventh';
    e = 4.0;
    ff = 7.0;
    g = 3.0;
    h = 2.0 / 3.0;
    x = a;
    [ x, ierror, k ] = t_family2 ( x, e, ff, g, h, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Eighth function in Traub, page 237.
%
  if ( 1 <= nder )

    name = 'Traub eighth';
    x = a;
    [ x, ierror, k ] = traub8 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Ninth function in Traub, page 237.
%
  if ( 1 <= nder )

    name = 'Traub ninth';
    x = a;
    [ x, ierror, k ] = traub9 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Type 1 family of iteration functions including
%  forms 10 and 11 of Traub, pages 237.
%
  if ( 1 <= nder )

    name = 'Traub type 1, form 10';
    rho = ( 1.0 - sqrt ( 5.0 ) ) / 2.0;
    e = 0.0;
    x = a;
    [ x, ierror, k ] = tt1f ( x, e, rho, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'Traub type 1, form 11';
    e = 1.0;
    x = a;
    [ x, ierror, k ] = tt1f ( x, e, rho, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Fourteenth function in Traub, page 237
%
  if ( 1 <= nder )

    name = 'Traub fourteenth';
    x = a;
    [ x, ierror, k ] = t14 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Fifteenth function in Traub, page 238.
%
  if ( 1 <= nder )

    name = 'Traub fifteenth';
    x = a;
    [ x, ierror, k ] = t15 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Sixteenth function in Traub, page 238.
%
  if ( 1 <= nder )

    name = 'Traub sixteenth';
    x = a;
    [ x, ierror, k ] = t16 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  A family of fourth order methods by Richard King.
%
  if ( 1 <= nder )

    name = 'King, BETA=0';
    beta = 0.0;
    x = a;
    [ x, ierror, k ] = king ( x, beta, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'King, BETA=1';
    beta = 1.0;
    x = a;
    [ x, ierror, k ] = king ( x, beta, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  if ( 1 <= nder )

    name = 'King, BETA=2';
    beta = 2.0;
    x = a;
    [ x, ierror, k ] = king ( x, beta, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Jarratt iterative technique.
%
  if ( 1 <= nder )

    name = 'Jarratt';
    x = a;
    [ x, ierror, k ] = jarratt ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Jarratt inverse-free iterative technique.
%
  if ( 1 <= nder )

    name = 'Jarratt inverse-free';
    x = a;
    [ x, ierror, k ] = jarratt2 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  4.  Multiple root methods, multiplicity given.\n' );
  fprintf ( 1, '\n' );
%
%  Script E2.
%
  if ( 1 <= nder )

    name = 'Script E2';
    x = a;
    [ x, ierror, k ] = script_e2 ( x, mult, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Script E3.
%
  if ( 2 <= nder )

    name = 'Script E3';
    x = a;
    [ x, ierror, k ] = script_e3 ( x, mult, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Script E4.
%
  if ( 3 <= nder )

    name = 'Script E4';
    x = a;
    [ x, ierror, k ] = script_e4 ( x, mult, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end
%
%  Star E - 1,1(f) function in Traub, page 235.
%
  name = 'Star E 1,1(f)';
  x = a;
  x1 = b;
  [ x, x1, ierror, k ] = te11f ( x, x1, mult, abserr, kmax, f );

  if ( ierror ~= 0 )
    fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
  else
    fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
  end
%
%  Routines which determine the root and multiplicity
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  5.  Multiple root methods,\n' );
  fprintf ( 1, '      multiplicity not given.\n' );
  fprintf ( 1, '\n' );
%
%  E2(U).
%
  if ( 2 <= nder )

    name = 'E2(U)';
    x = a;
    [ x, em, ierror, k ] = te2u ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d%15.6g\n', name, x, k, em );
    end

  end
%
%  Phi - 1,1(u) function in Traub, page 235.
%
  if ( 1 <= nder )

    name = 'Phi 1,1(U)';
    x = a;
    [ x, em, ierror, k ] = tphi1u ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d%15.6g\n', name, x, k, em );
    end

  end
%
%  Third function in Traub, page 235 bottom.
%
  if ( 1 <= nder )

    name = 'Traub third';
    x = a;
    [ x, em, ierror, k ] = tthip ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d%15.6g\n', name, x, k, em );
    end

  end
%
%  Van de Vel iteration.
%
  if ( 1 <= nder )
    name = 'Van de Vel';
    x = a;
    [ x, em, ierror, k ] = vandev ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d%15.6g\n', name, x, k, em );
    end

  end
%
%  Improved Van de Vel iteration.
%
  if ( 1 <= nder )
    name = 'Improved Van de Vel';
    x = a;
    [ x, em, ierror, k ] = vandev2 ( x, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d%15.6g\n', name, x, k, em );
    end

  end
%
%  Bisection methods
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  6. Bisection methods\n' );
  fprintf ( 1, '\n' );
%
%  Bisection methods may only be used if there is a change of sign.
%
  if ( ~ change )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'zoomin_all_test():\n' );
    fprintf ( 1, '  Bisection methods will not be called since\n' );
    fprintf ( 1, '  there is no change of sign interval.\n' );

  else
%
%  Bisection.
%
    name = 'Bisection';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = bisect ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end
%
%  Regula falsi.
%
    name = 'Regula falsi';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = regula ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end
%
%  Brent.
%
    name = 'Brent';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = brent ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end
%
%  Bisection + secant.
%
    name = 'Bisection + secant';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = rhein1 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end
%
%  Bisection + secant + inverse quadratic.
%
    name = 'Bisection + secant + inv quad';
    x = a;
    x1 = b;
    [ x, x1, ierror, k ] = rhein2 ( x, x1, abserr, kmax, f );

    if ( ierror ~= 0 )
      fprintf ( 1, '%30s%15.6g%8d error=%5d\n', name, x, k, ierror );
    else
      fprintf ( 1, '%30s%15.6g%8d\n', name, x, k );
    end

  end

  return
end
