new_raph( func, x0, maxiter, tol, data )

% This routine implements the Newton-Raphson's method of rootfinding.

% The advantage is that it can take any function for rootfinding
% The disadvantage is :
%    Currently it can deal only with single independent variable functions.
%    It does not deal with cases where the derivative of function becomes 0 in initial estimate itself.

% It can be directly invoked from MATLAB terminal or can be used as a sub-routine in some application.
% The following are the input parameters to the routine :
%    func : it is any function involving one independent variable of the form f( variables ) = 0.
%           i.e. all the terms of the function have to be transferred to LHS making RHS as 0.
%      x0 : initial guess of the root / starting value.
% maxiter : maximum number of iterations needed to be performed.
%     tol : allowable deviation / tolerance.
%    data : if the function 'func' has been defined in some '*.m' routine,
%           & is dependent on many other parameters,
%           those required data are packed into this input parameter.
%           in that case, the input function has to be appropriately built using data.
%           another option would have been to use varargin{:}; but that I leave for discretion of user.

function xc = new_raph( func, x0, maxiter, tol, data )

   h = 0.01;   % neighbourhood size
   iter = 0;   % iteration counter
   xc = x0;    % some random initialisations
   xp = 0;

   while( (abs(xp-xc) >= tol) && (iter <= maxiter) )
      xp = xc;

      % evaluate the function at given point
      fx = feval( func, xp, data );
      % evaluate the derivative using Central Difference Method
      fdashx = (feval( func, xp + h, data) + feval( func, xp - h, data))/ 2 / h;

      if( fdashx == 0 )
         iter = maxiter + 1;
         error( 'Derivative has become 0' );
      else
         xc = xp - fx / fdashx;
         iter = iter + 1;
      end

   end

   xc = ( xc + xp ) / 2;

return;