code: public static double cubic( double f0 , double g0 ,
  								double f1 , double alpha1 ,
  								double f2 , double alpha2 ) {
  		// Several different formulation were considered for solving this equation.
  		// See ExamineCubicInterpolateStability in benchmark directory
  		// Turns out that the straight forward implementation is just about the best.
  
  		double denominator = alpha1*alpha1*alpha2*alpha2*(alpha2-alpha1);
  		double a11 = alpha1*alpha1/denominator;
  		double a12 = -alpha2*alpha2/denominator;
  		double a21 = -alpha1*a11;
  		double a22 = -alpha2*a12;
  
  		double y1 = f2 - f0 - g0*alpha2;
  		double y2 = f1 - f0 - g0*alpha1;
  
  		double a = a11*y1 + a12*y2;
  		double b = a21*y1 + a22*y2;
  		
  	 	return (-b+Math.sqrt(b*b-3*a*g0))/(3.0*a);
  	}
docs: <p>
  Interpolates the next step using a cubic model.  Interpolation works by solving for 'a' and 'b' in
  the equation below. Designed to minimize the number of times the derivative
  needs to be computed.  Care has been taken reduce overflow/underflow by normalizing.
  </p>
  
  {@code
  &phi;(&alpha;) =  a*&alpha;<sup>3</sup>  + b*&alpha;<sup>2</sup> + &alpha;<sup>3</sup>  + &alpha;*&phi;'(0) + &phi;(0)
  }
  
  @param f0 Function value at f(0)
  @param g0 Derivative value at g(0)
  @param f1 Function value at f(a1)
  @param alpha1 value of a1
  @param f2 Function value at f(a2)
  @param alpha2 value of a2
  
  @return Interpolated step length