function x = bisection(fun,a,b)
fa = fun(a);
fb = fun(b);
 while (abs(fa-fb) > eps && (b-a) > eps)
    xx = (a+b)./2;
    fx = fun(xx);
    if (sign(fx) ~= sign(fa));
       b=xx;
       fb = fx;
    else
        a=xx;
        fa = fx;
    end
   end
x=xx;