function [x, k, vec] = mybis(f, h, a, b)
   %MYBISECT mybisect(f,a,b) reduces the interval [a,b] to a tiny
   % interval on which the function f(x) changes its sign.
   % Returns one end point of that interval, the number
   % of function evaluations, and the vector of approximations

   vec = [];
   fa = f(h, a);
   fb = f(h, b);
   if sign(fa) == sign(fb)
      error('f(x) must change sign on the interval [a,b].')
   end
   k = 2;
   
   while abs(a-b) > eps*abs(a)
       x = (a + b)/2;
       fx = f(h, x);
       k = k + 1;
       if sign(fx) == sign(fb)
           b = x;
           fb = fx;
       else
           a = x;
           fa = fx;
       end
       vec = [vec x];
	end
end

