function [f_p,x_p] = bisection(f,a,b)
Xa=a;
Xb=b;
Xm=(a+b)/2;
Fa=f(a);
Fb=f(b);
Fm=f(Xm);
err=abs(Xa-Xb);
while err > 1e-7
    Xl=(Xa+Xm)/2;
    Xr=(Xb+Xm)/2;
    Fl=f(Xl);
    Fr=f(Xr);
    F_max=max([Fa,Fb,Fm,Fl,Fr]);
    if F_max==Fa || F_max==Fl
        Xb=Xm;
        Xm=Xl;
        Fb=Fm;
        Fm=Fl;
    elseif F_max==Fm
        Xa=Xl;
        Xb=Xr;
        Fa=Fl;
        Fb=Fr;
    else
        Xa=Xm;
        Xm=Xr;
        Fa=Fm;
        Fm=Fr;
    end
    err=abs(Xa-Xb);
end
x_p=Xa;
f_p=f(x_p);

end