function retval = cubicmin(a,b,fa,fb,fpa,fpb)
% usage: x = cubicmin(a,b,fa,fb,fpa,fpb)
% description: returns the value x in the interval
% [a,b] at which the Hermite cubic interpolating 
% the boundary data f(a), f(b), f'(a), f'(b) 
% assumes a minimum value. We *assume* a<b.

% local variables:
% bma: b minus a
% d1,d2: temporaries for min formula
% absc: vector of candidate abscisssas
% ord: vector of candidate ordinates
% fab: divided difference f[a,b]
% faab: divided difference f[a,a,b]
% faabb: divided difference f[a,a,b,b]
% x: candidate abscissa of mininum

bma = b - a;
fab = (fb - fa)/bma;
d1 = fpa + fpb - 3*fab;
d2 = sqrt(d1*d1 - fpa*fpb);
absc = [a b];
ord = [fa,fb];
if (isreal(d2)) % any real roots?
  faab = (fab-fpa)/bma;
  faabb = (fpb-2*fab+fpa)/(bma*bma);
  if (abs(faabb) < 100*eps) % faabb nearly zero?
    if (abs(faab) < 100*eps) % faab nearly zero?
      x = a;
    else  
      x = a - fpa/(2*faab);
    end  
  else  
    x = b - bma*(fpb + d2 - d1)/(fpb - fpa + 2*d2);
  end
  x = max(min(x,b),a);
  absc = [a b x];
  ord = [fa,fb, fa + (x-a)*(fpa + (x-a)*(faab + (x-b)*faabb))];
end
retval = min(ord);
retval = find(ord == retval);
retval = absc(retval(1));