%=======================================================================
% ch11_bisection.m
%
% Find one root of f(x) = 0 by bisection on an initial interval.
%
% Origin:
%   Based on the bisection routine in Root5.m by Alexander Godunov
%   (February 2022). This standalone version uses the example,
%   stopping criterion, and status convention of ch11_bisection.f90.
%
% Example:
%   f(x) = sqrt(x) - cos(x) on [0,1].
%   Reference root = 0.64171437087288... .
%
% The program reports the iteration count, root estimate, and f(root).
% The requested tolerance applies to the width of the final bracket.
% A small bracket alone does not establish continuity or exclude a
% singularity, so inspect the residual and the function's domain.
%=======================================================================

clear; close all; clc;

% Initial bracket and requested maximum bracket width
x1  = 0.0;
x2  = 1.0;
eps = 1.0e-6;

fprintf('A single root of f(x)\n');
fprintf('Method - Bisection\n');

[root, flag] = bisection(@f, x1, x2, eps);

% Report a missing bracket or invalid input before printing a root
if flag == 0
    fprintf('invalid interval or no initial sign change\n');
    return
end

fprintf('iterations  = %d\n', abs(flag));
fprintf('tolerance   = %.5e\n', eps);
if flag < 0
    fprintf('iteration limit reached before the requested width\n');
end

fprintf('root        = %.5e\n', root);
% Check the residual of the reported estimate
fprintf('f(root)     = %.5e\n', f(root));


function [root, flag] = bisection(fun, x1, x2, eps)
%=======================================================================
% Locate one root of fun(x) on [x1,x2] by repeatedly halving the bracket.
%
% Input:
%   fun   - function handle, continuous on the initial bracket
%   x1    - left endpoint of the initial bracket
%   x2    - right endpoint of the initial bracket
%   eps   - requested maximum width of the final bracket (positive)
%
% Output:
%   root  - midpoint of the final bracket when flag is nonzero
%   flag  - convergence indicator and iteration count
%           flag > 0: width <= eps; flag is the iteration count
%           flag = 0: invalid interval/tolerance or no sign change
%           flag < 0: iteration limit reached; abs(flag) is the count
%
% Notes:
%   A sign change guarantees a root only when fun is continuous.
%   The maximum number of iterations is 200. If the requested width
%   is not reached by then, the estimate is not marked converged.
%   The algorithm uses bracket width, not |fun(root)|, to stop.
%=======================================================================

iter = 200;
flag = 0;
root = NaN;  % No estimate is returned when the initial bracket is invalid.

% Check that the inputs describe a usable initial bracket
if x1 >= x2 || eps <= 0
    return
end
if fun(x1)*fun(x2) > 0
    return
end

% Initialize the current bracket
a = x1;
b = x2;

% Repeatedly retain the half whose endpoints have opposite signs
for i = 1:iter
    % This midpoint form avoids overflow in a+b when both are large.
    c = a + (b-a)/2.0;
    if fun(a)*fun(c) <= 0
        b = c;
    else
        a = c;
    end

    % Stop when the current bracket is sufficiently narrow.
    if abs(b-a) <= eps
        break
    end
end

root = (b+a)/2.0;

% Report whether the requested bracket width was reached
if abs(b-a) <= eps
    flag = i;
else
    flag = -iter;
end
end


function y = f(x)
%-----------------------------------------------------------------------
% Example function. The chosen bracket must stay in the domain x >= 0.
% To try the alternative from the original file, use [-1,0].
%-----------------------------------------------------------------------
% y = x + cos(x);
y = sqrt(x) - cos(x);
end