function x2 = bisect(a,b,epsilon)
% bisection algorithm

x0=a; x1=b;
f0=f(x0); f1=f(x1);

% check initial bracket
if  f0*f1 >= 0 
   error('f0*f1 is non-negative!')
   end

% bisection loop
while abs(x1-x0) > epsilon
   x2 = (x0+x1)*0.5;
   f2 = f(x2);
disp([x0,x1,x2,f0,f1,f2]);
   if f2==0 
      return
      end
   if f0*f2<0    %root is to left of mid point
      x1=x2;
      f1=f2;
   else          %root is to right of mid point
      x0=x2;
      f0=f2;
   end
end

% final bisection to get x2 at centre of interval

x2=(x0+x1)*0.5;
