
function [ x, k, fout, xout ] = newton( func, dfunc, x0, atol, nmax )
% By Stathis Charalampidis
%
% On input:
%          a)  func: is the nonlinear function f(x).
%          b) dfunc: the derivative of f(x).
%          c)   x0 : initial guess.
%          d) atol : absolute tolerance.
%          e) nmax : maximum number of iterations allowed.
%
% On output:
%          a)    x : is the root of f(x).
%          b)    k : number of iterations required to achieve atol.
%          c)  fout: Residuals store at each k.

res = 10 * atol;  % Just to start while loop 
%xk - x0;                   % Compute the absolute error.
k = 0;                              % Iteration index.
iflag = 0;                         % For stopping purposes.
while abs(res)>atol            % Do-while loop.
    fxk = feval(func,x0);      % Compute the residual.
    fout(k+1) = fxk;            % Store the residual.
    xout(k+1) = x0;
    dfxk = feval(dfunc,x0);      % Evaluate the derivative of f(x).
    xk = x0 - fxk / dfxk;     % Newton step.
    res = xk - x0;             % Compute the "new" absolute error.
    x0 = xk;                   % Store it.
    k = k + 1;                 % Increase iteration index.
    if(k == nmax)              % Check whether we reacher the maximum
        iflag = -1;            % number of iterations!
        break;                 % If yes, then we stop Newton.
    end
end
if(iflag==0)                   % On successful exit, return the root!
    k = k + 1;                 % Increase iteration index.
    fout(k) = feval(func,xk);
    xout(k) = xk;
    x = xk;
else
    disp(' maximum number of iterations reached!');
    x = [];
    k = [];
end
end
