%=======================================================================
% ch11_secant.m
%
% Find one root of f(x) = 0 by the secant method.
%
% Origin:
%   MATLAB version of Secant.f90 by Alexander Godunov (January 2010).
%   The secant update and residual stopping test match ch11_secant.f90.
%   The two starting guesses are explicit. The old large-step midpoint
%   substitution is omitted because it was not a secant update.
%
% Example:
%   f(x) = x + cos(x), starting from x1 = -1 and x2 = -0.9.
%   Reference root = -0.73908513321516... .
%
% The program reports the iteration count, root estimate, and f(root).
% The secant method does not preserve a bracket. Starting values can
% affect convergence and may lead to a different root or failure.
%=======================================================================

clear; close all; clc;

% Two initial guesses and requested residual tolerance
x1  = -1.0;
x2  = -0.9;
eps = 1.0e-7;

fprintf('A single root of f(x)\n');
fprintf('Method - Secant\n');
fprintf('Initial guesses = %.3f, %.3f\n', x1, x2);

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

% Report invalid input or a zero secant slope before printing a root
if flag == 0
    fprintf('invalid tolerance/guesses or zero secant slope\n');
    return
end

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

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


function [root, flag] = secant(fun, x1, x2, eps)
%=======================================================================
% Locate one root using secant lines through successive iterates.
%
% Input:
%   fun   - function handle evaluating f(x)
%   x1    - first initial guess
%   x2    - second, distinct initial guess
%   eps   - requested maximum residual |fun(root)| (positive)
%
% Output:
%   root  - most recent secant estimate when flag is nonzero
%   flag  - convergence indicator and iteration count
%           flag > 0: |fun(root)| <= eps; flag is the iteration count
%           flag = 0: invalid input or zero secant slope
%           flag < 0: iteration limit reached; abs(flag) is the count
%
% Notes:
%   This is the ordinary secant method without damping or bracketing.
%   The maximum number of iterations is 200. An initial guess that
%   already satisfies the residual test is returned with flag=1 to
%   follow the positive-success convention, without an update.
%=======================================================================

iter = 200;
flag = 0;
root = NaN;  % No estimate is returned for invalid input or zero slope.
if eps <= 0 || x1 == x2
    return
end

a = x1;
b = x2;
fa = fun(a);
fb = fun(b);

% An accurate starting guess needs no secant update.
if abs(fa) <= eps
    root = a;
    flag = 1;
    return
end
if abs(fb) <= eps
    root = b;
    flag = 1;
    return
end

% Form the secant slope from the two most recent function values.
for i = 1:iter
    df = (fb-fa)/(b-a);
    if df == 0
        return
    end
    c = b - fb/df;
    fc = fun(c);

    % Stop on the residual, as in the original Fortran program.
    if abs(fc) <= eps
        root = c;
        flag = i;
        return
    end

    a = b;
    fa = fb;
    b = c;
    fb = fc;
    if a == b
        return
    end
end

root = c;
flag = -iter;
end


function y = f(x)
%-----------------------------------------------------------------------
% Example function. To use another problem, change this expression.
%-----------------------------------------------------------------------
y = x + cos(x);
end