# solving numerically a 2D integral by using simps and quad combined

I have the following function

$$f(x,y) =\left(\frac{1}{\exp(x-E_f)+1}-\frac{1}{\exp(x-y-E_f)+1} \right)\frac{1}{\sqrt{4-xy}}$$

I want to integrate it using Python 3. The domain of the $y$ variable is from sampled data hence I cannot used something like quad. The $x$ variable comes from analytical formula so I can generate whatever value I want for it. My question is how can use simps and quad combined to integrate this formula.