function r = bisection(f, a, b, tol)
	c = (a+b)/2;
%	hold on
	plot(c, f(c), 'o')
	if f(a)*f(b)>0 & abs(a-b)/2 < tol
		r = NaN;
		return
	elseif abs(a-b)/2 < tol
		r = c;
		return
	elseif (f(a)*f(c) < 0)
		r = bisection(f, a, c, tol);
	elseif (f(c)*f(b) < 0)
		r = bisection(f, c, b, tol);
	else
		r1 = bisection(f, a, c, tol);
		r2 = bisection(f, c, b, tol);
		if (r1 ~= NaN)
			r = r1
		elseif (r2 ~= NaN)
			r = r2
		else
			r = NaN;
		end
	end