def NR( f, f_, x0, TOL=1e-5, NMAX=100 ):
    """
    INPUT : - Function
            - Derivative of Function
            - Initial Guess
            - OPTIONAL: Tolerance
            - OPTIONAL: Maximum Number of Iterations

    RETURN: - Solution
            - Residual
            - Iterative Convergence
            - # of Iterations
    """
    n=1                                 # Counter
    res = []                            # Residual
    iter_conv = []                      # Iterative convergance

    while( n <= NMAX ):
        # Approximate value of x
        x1 = x0 - (f(x0)/f_(x0))
        res.append( f(x1) )             # Evaluate residual
        iter_conv.append( (x1-x0) )     # Evaluate iterative convergence
        if (abs(iter_conv[-1]) < TOL):
            # Return solution
            return (x1, res, iter_conv, n)
        
        else:
            x0 = x1                     # Update solution/guess
            n = n + 1                   # Increment counter

    return (False, res, iter_conv, n)