function [r, xk] = NRiter(f, fprime, x, tol, max_iter)
% NRiter: Newton-Raphson iterative method for root finding. Skeleton of
% a real implentation.
%
% [r, xk = NRiter(f, fprime, x0, tol)
%
% f        - the function
% fprime   - the differential of f
% x        - initial approximation
% tol      - tolerance
% max_iter - maximum number of iterations
%
% r        - the approximation of the root
% xk       - successive approximations of the root

if nargin < 4
    tol = 1e-6;
end
if nargin < 5
    max_iter = 100;
end

xk    = zeros(1, max_iter+1); 
xk(1) = x;

for k=1:max_iter
    % x = x - f(x) / fprime(x);
    xk(k+1) = xk(k) - f(xk(k)) / fprime(xk(k));
    if abs(f(xk(k+1))-f(xk(k))) < tol
        break;
    end
end

if k >= max_iter
    warning('Max. no. of iterations exhausted without reaching required tolerance.');
end

r = xk(k+1);
xk = xk(1:k+1);    % trim to the actual number of iterations performed
return

% Example:
%
%  f  = @(x) (x.^2 -4*sin(x));
%  fp = @(x) (2*x-4*cos(x));
%  NRiter(f, fp, 3)

