function [seq] = bisect(func, L, R, tol)
% implements the bisection method
% saves the sequence of midpoints as column vector

f=inline(func, 'x');

xl=L; % initial left boundary
xr=R; % initial right boundary

xm(1)=(xl+xr)/2;
k=1;

while (xr - xl)/2 > tol
    if f(xm(k)) == 0
        break % xm(k) is the root
    elseif sign( f(xm(k)) ) == sign( f(xl) )
        xl=xm(k); % move left boundary
    else 
        xr=xm(k); % move right boundary
    end
    k=k+1;
    xm(k)=(xl+xr)/2;
end

seq=xm.';

end

