%=======================================================================
% ch11_multiple_roots.m
%
% Find several distinct roots by scanning an interval and using
% bisection on each subinterval with a sign change.
%
% Origin:
%   MATLAB version of ch11_multiple_roots.f90 by Alexander Godunov.
%   Based on the default scan plus bisection method in BForce.f90.
%
% Example:
%   f(x) = exp(x) - x*sin(pi*x/2) on [-10.5,2].
%   The chapter example has six distinct roots in this interval.
%
% The tolerance applies to the final width of each root bracket.
% A sign-change scan can miss even-multiplicity roots, pairs of roots
% within one subinterval, and roots near a discontinuity. Inspect a
% plot and repeat with more subintervals when counting roots matters.
%=======================================================================

clear; close all; clc;

% Search interval, number of scan subintervals, and bracket-width tolerance
x1  = -10.5;
x2  = 2.0;
n   = 1000;
eps = 1.0e-7;

fprintf('Roots of f(x) by interval scanning and bisection\n');
fprintf('interval             = [%.3f, %.3f]\n', x1, x2);
fprintf('scan subintervals    = %d\n', n);
fprintf('bracket-width tol.   = %.5e\n', eps);

[roots, nfail] = multiple_roots(@f, x1, x2, n, eps);

if isempty(roots)
    fprintf('no roots detected by the scan\n');
else
    fprintf('number       root                  f(root)\n');
    for i = 1:numel(roots)
        fprintf('%5d   %18.10f   %16.8e\n', i, roots(i), f(roots(i)));
    end
end
if nfail > 0
    fprintf('brackets not refined to the requested width: %d\n', nfail);
end


function [roots, nfail] = multiple_roots(fun, x1, x2, n, eps)
%=======================================================================
% Scan [x1,x2] and bisect subintervals with opposite endpoint signs.
%
% Input:
%   fun     - function handle, continuous on the search interval
%   x1,x2   - endpoints of the search interval, with x1 < x2
%   n       - number of equal scan subintervals (positive integer)
%   eps     - maximum final bracket width (positive)
%
% Output:
%   roots   - detected roots, in increasing order
%   nfail   - number of sign-change brackets that reached the limit
%
% Notes:
%   A zero at a grid point is recorded once, including at x2.
%   A sign change only certifies a root when fun is continuous there.
%   The scan does not certify that all roots have been found.
%=======================================================================

roots = [];
nfail = 0;
if n <= 0 || n ~= floor(n) || x1 >= x2 || eps <= 0
    return
end

roots = zeros(1, n+1);  % At most one stored root per scan grid point.
nroots = 0;
dx = (x2-x1)/n;
a = x1;
fa = fun(a);

% Each grid point is checked once; a sign change brackets another root.
for j = 1:n
    b = x1 + j*dx;
    fb = fun(b);

    if fa == 0
        nroots = nroots+1;
        roots(nroots) = a;
    elseif fa*fb < 0
        [root, flag] = bisect_bracket(fun, a, b, eps);
        if flag > 0
            nroots = nroots+1;
            roots(nroots) = root;
        else
            nfail = nfail+1;
        end
    end

    a = b;
    fa = fb;
end

% The rightmost grid point is not a left endpoint in the loop.
if fa == 0
    nroots = nroots+1;
    roots(nroots) = a;
end
roots = roots(1:nroots);
end


function [root, flag] = bisect_bracket(fun, x1, x2, eps)
%=======================================================================
% Refine one sign-change bracket by bisection.
%
% Input:
%   fun     - function handle, continuous on [x1,x2]
%   x1,x2   - endpoints with opposite signs of fun
%   eps     - requested maximum final bracket width
%
% Output:
%   root    - midpoint of the final bracket
%   flag    - positive iteration count on convergence; -200 on limit
%
% Notes:
%   The stopping rule is bracket width, as in ch11_bisection.m.
%=======================================================================

iter = 200;
a = x1;
b = x2;
fa = fun(a);
flag = -iter;

for i = 1:iter
    c = a + (b-a)/2;
    fc = fun(c);
    if fa*fc <= 0
        b = c;
    else
        a = c;
        fa = fc;
    end
    if abs(b-a) <= eps
        flag = i;
        break
    end
end

root = a + (b-a)/2;
end


function y = f(x)
%-----------------------------------------------------------------------
% Example from Chapter 11. Replace this expression for another problem.
% Choose a scan interval on which the new function is continuous.
%-----------------------------------------------------------------------
pi_value = 3.141592653589793;
y = exp(x) - x*sin(pi_value*x/2);
end