function [ xi ] = bisection_min(f,I)
% function [ xi ] = bisection_min(f,I)
%  Demonstration of bisection minimisation method

nit=3;

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);

xx=linspace(x0,x4,101);
fx=f(xx);
Mf=max(fx);
mf=min(fx);
Df=Mf-mf;
ypl=mf;
plot(xx,fx)
axis([x0-0.1*d,x4+0.1*d,mf-0.75*Df,Mf+0.1*Df])
hold on
plot([x0,x4],(ypl-0.05*Df)*[1,1],'k')
plot([x0,x1,x2,x3,x4],ypl-0.05*Df,'r+')
text(x0-0.01*d,ypl-0.09*Df,'x_0')
text(x1-0.01*d,ypl-0.09*Df,'x_1')
text(x2-0.01*d,ypl-0.09*Df,'x_2')
text(x3-0.01*d,ypl-0.09*Df,'x_3')
text(x4-0.01*d,ypl-0.09*Df,'x_4')

for k=1:nit
    [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
    ypl=ypl-0.15*Df;
    plot([x0,x4],(ypl-0.05*Df)*[1,1],'k')
    plot([x0,x1,x2,x3,x4],ypl-0.05*Df,'r+')
    text(x0-0.01*d,ypl-0.09*Df,'x_0')
    text(x1-0.01*d,ypl-0.09*Df,'x_1')
    text(x2-0.01*d,ypl-0.09*Df,'x_2')
    text(x3-0.01*d,ypl-0.09*Df,'x_3')
    text(x4-0.01*d,ypl-0.09*Df,'x_4')
end
xi=(x3+x2)/2;

end

