function x = solve_zbrent(f, a, b, varargin)      % zero of f function in 1 dimension

tol  = 1e-14;

fa   = feval(f, a, varargin{:});
fb   = feval(f, b, varargin{:});


found = 0;
c     = b;
fc    = fb;
e     = abs(b-a);
d     = e;

while ~any(found)
    
   ifst = fb.*fc>0;
   c  = a.*ifst + c.*(~ifst);  %c is the new bracket
   fc = fa.*ifst + fc.*(~ifst);       
   e  = abs(b-a).*ifst + e.*(~ifst); %?
   d  = abs(b-a).*ifst + d.*(~ifst); %?
   
   ifst = abs(fc)>=abs(fb);
   a = a.*ifst + b.*(~ifst);
   b = b.*ifst+ c.*(~ifst);
   c = c.*ifst + a.*(~ifst);
   
   cond = abs(fc)>=abs(fb);
   fa = fa.*cond + fb.*~cond;   
   fb = fb.*cond + fc.*~cond;
   fc = fc.*cond + fa.*~cond;

   %check for convergence 
   xm   = 0.5*(c-b);                                                                                 % middle point            
   tol1 = 2*eps*abs(b) + 0.5*tol;
   found = min((abs(xm) < tol1)); % continue unless all are found

   %inverse quadratic interpolation
   invq = (abs(e) >= tol1).*(abs(fa) > abs(fb)); % First control to inverse quadratic, step has to be big enough     
   
   s   = fb./fa;                                                        
   r   = fb./fc;
   q   = fa./fc;
   
   p   = (a==c).*(2*xm.*s) + (a~=c).*(s.*((2*xm.*q).*(q-r) - (b-a).*(r-1)));
   q   = (a==c).*(1-s) + (a~=c).*((q-1).*(r-1).*(s-1));
   
   q   = q.*(p<=0) + (-1).*q.*(p>0); %check whether in bounds
   p   = abs(p);
    
   min1 = 3*xm.*q.*(q-r)-(b-a).*(r-1);
   min2 = abs(e.*q);
   invq  = invq & (2*p < min(min1,min2)) ; % if this fails, do bisection
   
   e_qi  = d;
   d_qi  = p./q;
   
   % bisection
   
   d_b = xm;
   e_b = d;
   
   % choose the step d and e of quadratic interpolation (if it passed the check, if not bisection)
   d_qi(isnan(d_qi)&(~invq))=0; %since NaN*0 = NaN and NaN +a = NaN, this makes d = NaN in the next step even when d_qi is irrelevant
   
   d = invq.*d_qi+(~invq).*d_b;
   e = invq.*e_qi+(~invq).*e_b;
            
   % this steps are independent of the step
   a  = b;
   fa = fb;
   b  = (b + d).*(abs(d)>tol1) + (b+sign(xm).*tol1).*(abs(d)<=tol1);
   fb = feval(f, b, varargin{:});
   
end


x = b;