code: public static double cubic2( double f0 , double g0 , double x0 ,
  								 double f1 , double g1 , double x1 )
  	{
  		double theta = 3.0*(f0-f1)/(x1-x0) + g0 + g1;
  		double s = Math.max(Math.abs(theta),Math.abs(g0));
  		s= Math.max(s,Math.abs(g1));
  		double gamma = s*Math.sqrt((theta/s)*(theta/s) - (g0/s)*(g1/s));
  		if( x1 < x0 )
  			gamma = -gamma;
  		double p = (gamma-g0) + theta;
  		double q = ((gamma-g0)+gamma) + g1;
  		
  		return x0 + (p/q)*(x1-x0);
  	}
docs: <p>
  Cubic interpolation using the function and derivative computed at two different points.  This particular
  implementation taken from [1] and appears to be designed to maximize stability.
  </p>
  <p>
  [1] MINPACK-2 source code http://ftp.mcs.anl.gov/pub/MINPACK-2/dcstep.f
  </p>
  
  @param f0 Value of f(x0)
  @param g0 Derivative f'(x0)
  @param x0 First sample point
  @param f1 Value of f(x1)
  @param g1 Derivative g'(x1)
  @param x1 Second sample point
  @return Interpolated point