import sympy as sp
integrand = psi[i]*psi[j]
I = sp.integrate(integrand,
(x, Omega[0][0], Omega[0][1]),
(y, Omega[1][0], Omega[1][1]))
# Fall back on numerical integration if symbolic integration
# was unsuccessful
if isinstance(I, sp.Integral):
integrand = sp.lambdify([x,y], integrand)
I = sp.mpmath.quad(integrand,
[Omega[0][0], Omega[0][1]],
[Omega[1][0], Omega[1][1]])