%
% Bisection for f(x) = 0
%
% f is contained in f.m
%
function [x]=bisec(a,b,tol)

format long;
%
% Plot the function for visualization
%
n = 100;
h = (b-a)/(n+1);
y = [a:h:b];
g = f(y);
figure(1)
hold;
plot(y,g,'r-');

% Bisection starts here 

x0 = a;
f0 = f(x0);

x1 = b;
f1 = f(x1);

plot(x0,f0,'bo')
plot(x1,f1,'bo')

while( abs(x1-x0) > tol ) 
    xm = 0.5*(x0+x1);
    fm = f(xm);
    plot(xm,fm,'bo')

    fprintf(1,'%13.7g\t %13.7g\t %13.7g\n', xm, fm, abs(x1-x0));
    pause;
    if( f0*fm <= 0 ) 
        x1 = xm;
        f1 = fm;
    else   
        x0 = xm;
        f0 = fm;
    end
end

x = xm;
hold;





