def bisection(f, a, b, max_iterations, error=0.0001):
    """Approximate solution of f(x)=0 on interval [a,b] by the bisection method.

    Parameters
    ----------
    f : function
        The function for which we are trying to approximate a solution f(x)=0.
    a,b : numbers
        The interval in which to search for a solution. The function returns
        None if f(a)*f(b) >= 0 since a solution is not guaranteed.
    max_iterations : (positive) integer
        The number of iterations to implement.
    error: (positive) float
        The max error that is accepted

    Returns
    -------
    x_N : number
        The midpoint of the Nth interval computed by the bisection method. The
        initial interval [a_0,b_0] is given by [a,b]. If f(m_n) == 0 for some
        midpoint m_n = (a_n + b_n)/2, then the function returns this solution.
        If all signs of values f(a_n), f(b_n) and f(m_n) are the same at any
        iteration, the bisection method fails and return None.

    Examples
    --------
    >>> f = lambda x: x**2 - x - 1
    >>> bisection(f,1,2,25)
    1.618033990263939
    >>> f = lambda x: (2*x - 1)*(x - 3)
    >>> bisection(f,0,1,10)
    0.5
    """
    if f(a) * f(b) >= 0:
        return None
    a_n = a
    b_n = b
    n = 0
    while error < abs(f(a_n) - f(b_n)) or n < max_iterations:
        m_n = (a_n + b_n) / 2
        f_m_n = f(m_n)
        if error > abs(f(a_n) - f(b_n)):
            return (a_n + b_n) / 2
        elif f(a_n) * f_m_n < 0:
            a_n = a_n
            b_n = m_n
        elif f(b_n) * f_m_n < 0:
            a_n = m_n
            b_n = b_n
        elif f_m_n == 0:
            return m_n
        else:
            return None
    return (a_n + b_n) / 2
