function p=bisection(f,a,b,tol,max_iterations)

% Bisection: See section 2.1.
for i=1:max_iterations
    p=(a+b)/2;
    if p-a<tol
        return
    end
    if f(a)*f(p)>0
        a=p;
    else
        b=p;
    end
end
error('Exceeded maximum number of iterations');