simpson_cosmo

integral= simpson(function, a, b)
or integral= simpson(function, a, b, epsilon)
returns the integral of FUNCTION(x) from A to B. If EPSILON is
given, Simpson's rule is refined until that fractional accuracy
is obtained. EPSILON defaults to 1.e6.
If the notvector= keyword is supplied and nonzero, then FUNCTION
may not be called with a list of x values to return a list of
results. By default, FUNCTION is assumed to be a vector function.
If the function is very smooth, romberg may work better.
interpreted function, defined at ./cosmo.i line 331
