function zero_muller_test ( )

%*****************************************************************************80
%
%% zero_muller_test() tests zero_muller().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../zero_muller' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'zero_muller_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version );
  fprintf ( 1, '  Test zero_muller(), which uses Muller''s method,\n' );
  fprintf ( 1, '  with complex arithmeic, to solve a nonlinear equation.\n' );

  test01 ( );
  test02 ( );
  test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'zero_muller_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  rmpath ( '../zero_muller' );

  return
end
function test01 ( )

%*****************************************************************************80
%
%% test01() tests zero_muller() on F(X) = X*X+9.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test01():\n' );
  fprintf ( 1, '  Demonstrate zero_muller() on F(X) = X*X+9.\n' );

  fatol = 1.0E-05;
  itmax = 10;
  x1 = complex ( 1.0, 0.0 );
  x2 = complex ( 0.0, 1.0 );
  x3 = complex ( 0.5, 0.5 );
  xatol = 1.0E-05;
  xrtol = 1.0E-05;

  [ xnew, fxnew ] = zero_muller ( @func01, fatol, itmax, x1, x2, x3, ...
    xatol, xrtol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     X   = %20.10f,%20.10f\n', real ( xnew ), imag ( xnew ) );
  fprintf ( 1, '  with function value F(X):\n' );
  fprintf ( 1, '     FX   = %20.10f,%20.10f\n', real ( fxnew ), imag ( fxnew ) );
  fprintf ( 1, '   ||FX|| = %20.10f\n', abs ( fxnew ) );

  return
end
function test02 ( )

%*****************************************************************************80
%
%% test02() tests zero_muller() on F(X) = (X*X+4) * (X-10) * (X+20)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test02():\n' );
  fprintf ( 1, '  Demonstrate zero_muller() on F(X) = (X*X+4)*(X-10)*(X+20).\n' );

  fatol = 1.0E-05;
  itmax = 10;
  x1 = complex ( 1.0, 0.0 );
  x2 = complex ( 0.0, 1.0 );
  x3 = complex ( 0.5, 0.5 );
  xatol = 1.0E-05;
  xrtol = 1.0E-05;

  [ xnew, fxnew ] = zero_muller ( @func02, fatol, itmax, x1, x2, x3, ...
    xatol, xrtol );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     X   = %20.10f,%20.10f\n', real ( xnew ), imag ( xnew ) );
  fprintf ( 1, '  with function value F(X):\n' );
  fprintf ( 1, '     FX   = %20.10f,%20.10f\n', real ( fxnew ), imag ( fxnew ) );
  fprintf ( 1, '   ||FX|| = %20.10f\n', abs ( fxnew ) );

  return
end
function test03 ( )

%*****************************************************************************80
%
%% test03() tests zero_muller() on Zhelyazkov's function
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test03():\n' );
  fprintf ( 1, '  Demonstrate zero_muller() on Zhelyazkov''s function.\n' );

  fatol = 1.0E-07;
  itmax = 10;

  for i = 1 : 2
%
%  First set of starting points.
%  Result is X = ( 1.5705798926, 0.0 )
%
    if ( i == 1 )
      x1 = complex ( 1.0, 0.0 );
      x2 = complex ( 0.0, 1.0 );
      x3 = complex ( 0.5, 0.5 );
%
%  Second set of starting points.
%  Result is X = ( -0.5802520567, 0.0 ).
%
    elseif ( i == 2 )
      x1 = complex (  0.0, 1.0 );
      x2 = complex (  1.0, 2.0 );
      x3 = complex ( -1.0, 2.0 );
    end

    xatol = 1.0E-07;
    xrtol = 1.0E-07;

    [ xnew, fxnew ] = zero_muller ( @func03, fatol, itmax, x1, x2, x3, ...
      xatol, xrtol );

    fprintf ( 1, '\n' );
    fprintf ( 1, '     X   = %20.10f,%20.10f\n', real ( xnew ), imag ( xnew ) );
    fprintf ( 1, '  with function value F(X):\n' );
    fprintf ( 1, '     FX   = %20.10f,%20.10f\n', real ( fxnew ), imag ( fxnew ) );
    fprintf ( 1, '   ||FX|| = %20.10f\n', abs ( fxnew ) );

  end

  return
end
function fx = func01 ( x )

%*****************************************************************************80
%
%% func01() evaluates F(X) = X*X+9.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    complex X: the evaluation point.
%
%  Output:
%
%    complex FX: the function value.
%
  fx = x * x + 9.0;
 
  return
end
function fx = func02 ( x )

%*****************************************************************************80
%
%% func02() evaluates F(X) = (X*X+4)*(X-1)*(X+2).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    complex X: the evaluation point.
%
%  Output:
%
%    complex FX, the function value.
%
  fx = ( x * x + 4.0 ) * ( x - 10.0 ) * ( x + 20.0 );
 
  return
end
function fz = func03 ( z )

%*****************************************************************************80
%
%% func03() evaluates Zhelyazkov's function.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    complex Z: the evaluation point.
%
%  Output:
%
%    complex FZ: the function value.
%
  epsilon = complex ( 0.4, 0.0 );
  eta = complex ( 0.64, 0.0 );
  me = 0.384;
  mo = 0.5;
  one = complex ( 1.0, 0.0 );
  x = 0.5;

  ok = z - me / sqrt ( eta );
  of = z - mo;

  a = of * of + ( ok * ok ) * eta * tanh ( x );

  b = ( of - ok * eta ) / ( of - ok * eta * eta );

  fz = of * of - one + ( eta * ok * ok - one ) * ...
    tanh ( x ) - x * x * epsilon * epsilon * a * b;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the 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

