# Estimate the derivatives at three points {x01,x12,x23}:
      x01 = (x0 + x1)/2; d01 = (y1 - y0)/(x1 - x0) 
      x12 = (x1 + x2)/2; d12 = (y2 - y1)/(x2 - x1) 
      x23 = (x2 + x3)/2; d23 = (y3 - y2)/(x3 - x2) 
      # Estimate the second derivatives at two points {x012,x123}:
      x012 = (x01 + x12)/2; dd012 = (d12 - d01)/(x12 - x01) 
      x123 = (x12 + x23)/2; dd123 = (d23 - d12)/(x23 - x12) 
      # Interpolate the second derivative at {x}:
      s = (x - x012)/(x123 - x012)
      dd = (1 - s)*dd012 + s*dd123
      # dd_cor = -A*A*math.sin(A*x)
      # sys.stderr.write("dd = %.8f dd_cor = %.8f\n" % (dd, dd_cor))
      # Apply the quadratic correction:
      dy2 = -dd*r*(1-r)*(x2-x1)**2/2
      y += dy2