function x0 = bisectionRoot_JMS( y, x1, x2, dxMax )

% Find a root x0 of function func, bracketed between x1 and x2
% Input:
%   y(x)   : function whose root is to be found
%   x1, x2 : lower and upper bounds of root
%   dxMax  : max. error in value of root x0 (optional)
% Output:
%   x0     : root of y(x), i.e. y(x0)=0
% Usage:
%   x0 = bisectionRoot( @myFunc, 0, 2, 1.e-8 )
% J.M.Soler. Oct.2011

% Set dxMax, if not present
if ~exist('dxMax')
    dxMax = 1.e-12;   % default value
end

% Check bounds
y1 = y(x1);
y2 = y(x2);
if (y1*y2>0)
    error('bisectionRoot ERROR: x1,x2 not bracket of root')
end

% Bisection iteration
while (abs(x2-x1)>dxMax)
    x0 = (x1+x2)/2;
    y0 = y(x0);
    if (y0*y1>0) 
        x1 = x0;
        y1 = y0;
    else
        x2 = x0;
        y2 = y0;
    end
end % while
x0 = (x1*y2-x2*y1) / (y2-y1);

end % function
