function root = newtonraphson( func, dfunc, xr, es, maxit)

% NEWTONRAPHSON: root = newtraph(func,dfunc,xguess,es,maxit):
%   uses Newton-Raphson method to find the root of a function
% input:
%   func   = name of function 
%   dfunc  = name of derivative of function 
%   xguess = initial guess
%   es     = (optional) stopping criterion (%)
%   maxit  = (optional) maximum allowable iterations
% output:
%   root = real root

% if necessary, assign default values

if nargin<5, maxit=50; end     %if maxit blank set to 50
if nargin<4, es=0.001; end     %if es blank set to 0.001

% Newton-Raphson Iterates

iter = 0;
while (1)
  xrold = xr;
  xr = xr - func(xr)/dfunc(xr);
  iter = iter + 1;
  if xr ~= 0, ea = abs((xr - xrold)/xr) * 100; end
  if ea <= es | iter >= maxit, break, end
end
root = xr;