# How to implement Simpson's rule for double integral (without numeric limits of first integral)

I want to use Simpson's rule to evaluate the following double integral: $$\int_{a}^{b}\left|\int_{0}^{z}x\cdot \mathrm{erf}(x-10)\cdot J_{0}(x) \mathrm{dx}\right|^{2}\exp(-0.5*(z-40)^2)\mathrm{dz}$$

Usually, if both the limits were numeric I would get away with just applying simps from scipy.integrate twice, but here I have no idea how to proceed.

Your integral can be written as $$\int_a^b f(z)^2 h(z) dz \approx \sum_q f(z_q)^2 h(z_q) w_q$$ where I wrote some quadrature approximation. Now do another quadrature for $$f(z)$$ $$f(z) = \int_0^z g(x) dx \approx \sum_q g(x_q) v_q$$ where $$x_q \in [0,z]$$ are some quadrature nodes and $$v_q$$ are corresponding weights.
Write a function to compute $$f(z)$$ that does the second quadrature. Then use that do the first quadrature.