Problems integrating using quadpy

I am trying to use quadpy as I want to do 2D numerical vector integration for 2D integrals. To see how fast quadpy works, I wanted to test it and compare it with scipy 1D vector integration. Thus I wrote a simple code:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(q * t)


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y)
print(y1)

Having no experience with quadpy, I get the following errors:

Traceback (most recent call last):
  File "test.py", line 8, in <module>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 44, in quad
    g, [a, b], eps_abs=epsabs, eps_rel=epsrel, max_num_subintervals=limit
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_tools.py", line 42, in integrate_adaptive
    range_shape=range_shape,
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_gauss_kronrod.py", line 124, in _gauss_kronrod_integrate
    fx_gk = numpy.asarray(f(sp))
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 41, in g
    return f(x, *args)
  File "test.py", line 8, in <lambda>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
ValueError: operands could not be broadcast together with shapes (500,) (15,)

You're missing one multiply.outer:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(np.multiply.outer(q, t))


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y - y1)