function [ a, b, mx, my, nf ] = optisection ( f, a, b, n, x_tol, c )

%*****************************************************************************80
%
%% optisection() seeks a minimizer of a unimodal function in [a,b].
%
%  Discussion:
%
%    A strict unimodal function f(x) in [a,b] has a minimizer c such that
%    f(x) strictly decreases between a and c, and strictly increases between
%    c and b.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 May 2026
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Vladimir Kodnyanko,
%    Description and software implementation of optisection method
%    for minimizing unimodal functions.
%
%    Vladimir Kodnyanko,
%    Optimal section method for minimizing unimodal functions,
%    Computational Mathematics and Modeling,
%    Volume 36, pages 235-252, 16 October 2025.
%
%  Input:
%
%    function handle f: the function to be minimized.
%
%    real a, b: the left and right endpoints.
%
%    integer n: the maximum number of iterations allowed.
%    n must be at least 1.
%
%    real x_tol: a tolerance for the width of the x interval.
%
%    real c: the proportionality coefficient.  A suggested value is 0.2.
%
%  Output:
%
%    real a, b: the current interval containing the minimizer.
%
%    real mx, my: the estimated minimizer, and the minimizing value.
%
%    integer nf: the number of function evaluations.
%
  nf = 0;
%
%  At the end of each iteration, 
%  the minimizer is estimated to occur at a <= mx <= b, with value f(mx) = my.
%
  for it = 1 : n
%
%  On step 0, the interval is [a,b].
%
    if ( it == 1 )

      mx = ( a + b ) / 2.0;
      my = f ( mx );
%
%  If [ a, mx ] is longer than [ mx, b ]:
%
    elseif ( b - mx <= mx - a )

      px = c * a + ( 1.0 - c ) * mx;
      py = f ( px );

      if ( my < py )
        a = px;
      else
        b = mx;
        mx = px;
        my = py;
      end
%
%  [ mx, b ] is longer than [ a, mx ].
%
    else

      px = c * b + ( 1.0 - c ) * mx;
      py = f ( px );

      if ( my < py )
        b = px;
      else
        a = mx;
        mx = px;
        my = py;
      end

    end
%
%  The function has been evaluated one more time.
%
    nf = nf + 1;
%
%  Decide if we can stop.
%
    if ( b - a <= x_tol )
      return
    end

  end

  return
end

