I have the following integral that I would like to evaluate:

`A[x_, y_] := 2/\[Pi] ((-1)^x y^x)/(2 - y)^(1 + x); Pintegrand[m_, n_, r1_, r2_, \[Theta]1_, \[Theta]2_] := (A[m, \[Tau]] A[n, \[Sigma]])/(\[Pi]^2) r1 r2 LaguerreL[m, (4 r1^2)/(2 - \[Tau])] LaguerreL[n, (4 r2^2)/(2 - \[Sigma])] Exp[-((2 \[Tau] r1^2)/(2 - \[Tau]))] Exp[-((2 \[Sigma] r2^2)/(2 - \[Sigma]))] Exp[- (1/2) (-R2 + r2 Cos[\[Theta]2])^2 + (-I1 + r1 Sin[\[Theta]1]) (-((b (-R1 + r1 Cos[\[Theta]1]))/(-b^2 + a c)) + (a (-I1 + r1 Sin[\[Theta]1]))/(-b^2 + a c)) + (-R1 + r1 Cos[\[Theta]1]) ((c (-R1 + r1 Cos[\[Theta]1]))/(-b^2 + a c) - (b (-I1 + r1 Sin[\[Theta]1]))/(-b^2 + a c)) + (-I2 + r2 Sin[\[Theta]2])^2] // Simplify; P11 = Integrate[Pintegrand[1, 1, r1, r2, \[Theta]1, \[Theta]2], {r1, 0, \[Infinity]}, {\[Theta]1, 0, 2 \[Pi]}, {r2, 0, \[Infinity]}, {\[Theta]2, 0, 2 \[Pi]}] `

The integrand function is the integrand of the integral, which is a Gaussian function weighted with the Laguerre polynomials (which is what is making the calculation lengthy).

I would like to generate a list of results with different values of $ m$ and $ n$ – that is P01, P10, P11, P12, P21, P22, P13, etc (first number corresponds to value of $ m$ and second to $ n$ ), with both $ m$ and $ n$ integers that each go from 0 to 10. I would like to use the result to create a contour plot of the result.

There are two post that request a similar objective, such as this and this. However the solutions appear to be specific to the integral in question and make use of different strategies, such as `Parallelize`

, `ParallelTable`

, and `ParallelMap`

.

How can I most efficiently parallelise the numerical integral to generate a contour plot? Any help is appreciated!