function root = find_root_newtonraphson (fun, deriv, xinit, desired_acc)

%NEWTONRAPHSON     Determines the root of a function
%   ROOT = NEWTONRAPHSON (FUN, DERIV, XINIT, DESIRED_ACC) determines the
%   ROOT of the function defined by the function handle FUN with the
%   corresponding derivative defined by the function handle DERIV.  XINIT
%   is use as an initial guess.  The function iterates until the desired
%   accuracy DESIRED_ACC has been reached.  An error is returned when a
%   root cannot be determined.
%
%   Inputs:
%       FUN:    a function handle to the function for which the root is
%               solved
%       DERIV:  a function handle to the derivative of FUN
%       XINIT:  the initial guess for the root
%       DESIRED_ACC:  an upper bound on the approximate absolute error
%
%   Outputs:
%       ROOT:   an approximation of the root within the desired accuracy
%
%   Assumptions:
%       - if the function does not converge after 100 iterations, it is
%       assumed that a root cannot be determined
%       - DERIV represents an accurate implementation of the derivative of
%       FUN
%
%   Example Usage:
%       >> root = newtonraphson(@sin, @cos, 1, 1e-6);
%
max_iter = 100;
root_previous = xinit;
for i=1:max_iter
    f_val = feval(fun,root_previous);
    deriv_val = feval(deriv,root_previous);
    if deriv_val == 0
        error('did not converge -- zero derivative');
    end
    root = root_previous - f_val/deriv_val;
    if abs(root - root_previous) < desired_acc
        return
    end
    root_previous = root;
end
error('did not converge -- max iterations exceeded');