%=======================================================================
% ch11_newton.m
%
% Find one root of f(x) = 0 by Newton's method.
%
% Origin:
%   MATLAB version of Newton1.f90 by Alexander Godunov (January 2010).
%   The Newton update and residual stopping test match ch11_newton.f90.
%   The old large-step substitution is omitted because it was neither
%   a Newton step nor a bisection step.
%
% Example:
%   f(x) = x + cos(x), starting from x0 = -1.
%   Reference root = -0.73908513321516... .
%
% The program reports the iteration count, root estimate, and f(root).
% Newton's method is local: a different initial guess may converge to
% another root or fail. Interpret a small residual in the context of
% the function and the physical problem.
%=======================================================================

clear; close all; clc;

% Initial guess and requested residual tolerance
x0  = -1.0;
eps = 1.0e-7;

fprintf('A single root of f(x)\n');
fprintf('Method - Newton\n');
fprintf('Initial guess = %.3f\n', x0);

[root, flag] = newton1(@f, x0, eps);

% Report invalid input or a zero derivative before printing a root
if flag == 0
    fprintf('invalid tolerance or zero derivative at an iterate\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, 0));


function [root, flag] = newton1(fun, x0, eps)
%=======================================================================
% Locate one root near x0 using Newton's update x_new=x-f(x)/f'(x).
%
% Input:
%   fun   - function handle evaluating f(x) for k=0 and f'(x) for k=1
%   x0    - initial guess
%   eps   - requested maximum residual |f(root)| (positive)
%
% Output:
%   root  - most recent Newton estimate when flag is nonzero
%   flag  - convergence indicator and iteration count
%           flag > 0: |f(root)| <= eps; flag is the iteration count
%           flag = 0: invalid tolerance or zero derivative
%           flag < 0: iteration limit reached; abs(flag) is the count
%
% Notes:
%   This is the ordinary Newton method without damping or bracketing.
%   The maximum number of iterations is 200. If the initial guess
%   already satisfies the residual test, flag=1 denotes success even
%   though no Newton update was needed.
%=======================================================================

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

% A sufficiently accurate initial guess needs no Newton update.
fx = fun(x, 0);
if abs(fx) <= eps
    root = x;
    flag = 1;
    return
end

% Repeatedly use the value and derivative at the current estimate.
for i = 1:iter
    dfx = fun(x, 1);
    if dfx == 0
        return
    end

    xnew = x - fx/dfx;
    fx = fun(xnew, 0);

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

root = xnew;
flag = -iter;
end


function y = f(x, k)
%-----------------------------------------------------------------------
% Evaluate the example function for k=0 and its derivative for k=1.
% To use another function, change both expressions consistently.
%-----------------------------------------------------------------------
if k == 0
    y = x + cos(x);
elseif k == 1
    y = 1 - sin(x);
else
    error('The function selector k must be 0 or 1.');
end
end