ble extrapolation is not used (QUADPACK uses it
for infinite intervals). This is because the algorithm here is
supposed to work on vector-valued functions, in an user-specified
norm, and the extension of the epsilon algorithm to this case does
not appear to be widely agreed. For max-norm, using elementwise
Wynn epsilon could be possible, but we do not do this here with
the hope that the epsilon extrapolation is mainly useful in
special cases.

References
----------
[1] R. Piessens, E. de Doncker, QUADPACK (1983).

Examples
--------
We can compute integrations of a vector-valued function:

>>> from scipy.integrate import quad_vec
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> alpha = np.linspace(0.0, 2.0, num=30)
>>> f = lambda x: x**alpha
>>> x0, x1 = 0, 2
>>> y, err = quad_vec(f, x0, x1)
>>> plt.plot(alpha, y)
>>> plt.xlabel(r"$\alpha$")
>>> plt.ylabel(r"$\int_{0}^{2} x^\alpha dx$")
>>> plt.show()

When using the argument `workers`, one should ensure
that the main module is import-safe, for instance
by rewriting the example above as:

.. code-block:: python

    from scipy.integrate import quad_vec
    import numpy as np
    import matplotlib.pyplot as plt

    alpha = np.linspace(0.0, 2.0, num=30)
    x0, x1 = 0, 2
    def f(x):
        return x**alpha

    if __name__ == "__main__":
        y, err = quad_vec(f, x0, x1, workers=2)
N