function [ xi, nit ] = bisecmin(f,I,tol)
% function [ xi ] = bisecmin(f,I,tol)
%  Bisection minimisation method

if nargin==2
    tol=0.000001;
end

x0=I(1);
x4=I(2);
d=(x4-x0);
x2=(x0+x4)/2;
x1=(x0+x2)/2;
x3=(x2+x4)/2;
f1=f(x1);f2=f(x2);f3=f(x3);

nit=1;

while abs(x3-x1)>tol
    [p,mk]=min([f1,f2,f3]);
    if mk==1
        x4=x2;
        x2=x1; f2=f1;
        x1=(x0+x2)/2;
        x3=(x2+x4)/2;
        f1=f(x1);f3=f(x3);
    elseif mk==2
        x0=x1;x4=x3;
        x1=(x0+x2)/2;
        x3=(x2+x4)/2;
        f1=f(x1);f3=f(x3);
    else
        x0=x2;
        x2=x3;f2=f3;
        x1=(x0+x2)/2;
        x3=(x2+x4)/2;
        f1=f(x1);f3=f(x3);
    end
    nit=nit+1;
end
xi=(x3+x2)/2;

end

