f must return the value of the integrand at a given point.
Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by
d01pac. If your code inadvertently
does return any NaNs or infinities,
d01pac is likely to produce unexpected results.