function z=bisection(f,lb,ub,eps)
%INPUT
%================
%f ................... a scalar function
%lb ................. the initial lower bound
%ub ................ the initial upper bound
%eps .............. tolerance parameter
%OUTPUT
%================
% z ................. a root of the equation f(x)=0

if (f(lb)*f(ub)>0)
    error('f(lb)*f(ub)>0')
end

iter=0;
while (ub-lb>eps)
    z=(lb+ub)/2;
    iter=iter+1;
    if(f(lb)*f(z)>0)
        lb=z;
    else
        ub=z;
    end
    fprintf('iter_numer = %3d current_sol = %2.6f \n',iter,z);
end


