function range=findbracket(f,x0,dx,Nmax)
%finds interval [a,b] around x0 s.t. f(a)*f(b)<0
%steps: a=b=x0, dx=0.001 -> a=a-dx, check fa*fb<0 -> b=b+dx, check again ->
%dx=2*dx and reiterate from step 2.

a=x0;
b=x0;
if nargin == 2
    dx=0.001;
    Nmax = 10^6;
end
if nargin == 3
   Nmax = 10^6; 
end

iter = 1;
while iter <= Nmax
    a = a-dx;
    if f(a)*f(b) < 0
        break, end;
    b = b+dx;
    if f(a)*f(b) < 0
        break, end;
    dx = 2*dx; 
    iter = iter + 1;
end
range = [a,b];
if iter > Nmax
   disp('May not have found an appropriate interval...') 
end
end

