title: What would be the computationally faster way to implement this 2D numerical integration?
body: <p>I am interested in doing a 2D numerical integration. Right now I am using the <code>scipy.integrate.dblquad</code> but it is very slow. Please see the code below. My need is to evaluate this integral 100s of times with completely different parameters. Hence I want to make the processing as fast and efficient as possible. The code is:</p>
  
  <pre><code>import numpy as np
  from scipy import integrate
  from scipy.special import erf
  from scipy.special import j0
  import time
  
  q = np.linspace(0.03, 1.0, 1000)
  
  start = time.time()
  
  def f(q, z, t):
      return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
          -0.5 * ((z - 40) / 2) ** 2)
  
  
  y = np.empty([len(q)])
  for n in range(len(q)):
      y[n] = integrate.dblquad(lambda t, z: f(q[n], z, t), 0, 50, lambda z: 10, lambda z: 60)[0]
  
  end = time.time()
  print(end - start)
  
  </code></pre>
  
  <p>Time taken is </p>
  
  <pre><code>212.96751403808594
  </code></pre>
  
  <p>This is too much. Please suggest a better way to achieve what I want to do. I tried to do some search before coming here, but didn't find any solution. I have read <code>quadpy</code> can do this job better and very faster but I have no idea how to implement the same. Please help.  </p>