Integrate function using adaptive Simpson quadrature.
pcoriaccel_integrate_simp(integrand, lo, hi, tol = 1.490116e-08)integration result, list with elements $Q (the integral estimate), $fcnt (the number
of function evaluations), and $estim.prec (a (pessimistic) estimate of the precision).
The integrand, must take scalar argument, may return scalar, vector, or matrix.
Lower integration bound
Upper integration bound
Tolerance for integration, default .Machine$double.eps^(1/2)