Lets say you want to evaluate the double integral $$\int_1^3 \int_0^2 (y^3 - xy) dy dx.$$ This can be done using scipy's "dblquad" function, as is shown below.
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi
def f(x,y):
return y**3 - x*y
lim = 5
xv, yv = lim*np.linspace(-1,1,500), lim*np.linspace(-1,1,500)
X,Y = np.meshgrid(xv,yv)
plt.figure(dpi=150)
plt.imshow(f(X,Y), extent=(-lim,lim,-lim,lim))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.plot(np.linspace(1,3,100), 0.0*np.ones(100), 'r-', lw=3)
plt.plot(np.linspace(1,3,100), 2.0*np.ones(100), 'r-', lw=3)
plt.plot(np.ones(100), np.linspace(0,2,100), 'r-', lw=3)
plt.plot(3.0*np.ones(100), np.linspace(0,2,100), 'r-', lw=3)
plt.show()
# You need to be careful with the order of the limits to get the correct answer
I = spi.dblquad(f, 0,2, lambda x: 1, lambda x: 3)
print(I)
(2.220446049250313e-16, 8.785604931891442e-14)
We know, from evaluating this analytically that the answer is zero.