Question: I have a function that takes as an input a scalar value, and outputs a matrix (following is an MRE with size 2x2, my actual function has problem specific matrix dimensions, like say 100x100). I'd like to integrate said function over a range. import numpy as np from scipy.integrate import quad def func(x): return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]]) quad(func, 0, np.pi)  This gives me the following error : Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/scipy/integrate/_quadpack_py.py", line 465, in quad retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/scipy/integrate/_quadpack_py.py", line 577, in _quad return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: only size-1 arrays can be converted to Python scalars  I wanted to avoid looping over each individual function, so I vectorized the integration operation, in the hopes that it would accommodate this sort of a function : grid_quad = np.vectorize(quad) grid_quad(func, 0, np.pi)  I still get the same error : Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/numpy/lib/function_base.py", line 2329, in __call__ return self._vectorize_call(func=func, args=vargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/numpy/lib/function_base.py", line 2407, in _vectorize_call ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/numpy/lib/function_base.py", line 2367, in _get_ufunc_and_otypes outputs = func(*inputs) ^^^^^^^^^^^^^ File "/scipy/integrate/_quadpack_py.py", line 465, in quad retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/scipy/integrate/_quadpack_py.py", line 577, in _quad return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: only size-1 arrays can be converted to Python scalars  When I explicitly wrote out my integrand as an array of arrays, this method worked efunc = [[np.cos, np.sin],[np.sin, np.cos]] grid_quad(efunc, 0, np.pi)  and gave the expected output (array([[4.92255263e-17, 2.00000000e+00], [2.00000000e+00, 4.92255263e-17]]), array([[2.21022394e-14, 2.22044605e-14], [2.22044605e-14, 2.21022394e-14]]))  which is a matrix that contains the values obtained by integrating the "function" associated with each matrix position. In pen and paper terms, I have a matrix whose entries are functions that depend on a variable. I would like to integrate each one of these functions over an interval, and receive as output a matrix of the definite integral values. Is there a "(num)pythonic" way to achieve this?

Answer: Try scipy.integrate.quad_vec. import numpy as np from scipy.integrate import quad_vec def func(x): return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]]) quad_vec(func, 0, np.pi) # (array([[2.22044605e-16, 2.00000000e+00], # [2.00000000e+00, 2.22044605e-16]]), # 1.3312465980656846e-13)  The first output is your result, the second is the error estimate.