How to solve this complicated double integral problem? I have written this in way to make it as much as possible non-confusing. I will start describing my problem and I will walk you through my question,
I have a double integration which I am trying to solve of the following form, 
$$F=\int\limits_{0<\gamma_1<\gamma_2<+\infty} \mathcal{L}\bigl(\sum_{i=1}^2\gamma_i^{-1} \bigl)\ g(\gamma_1,\gamma_2)\ d\gamma_1 d\gamma_2$$
where  $$\mathcal{L}(t) = \text{exp}\left(- \int_{\gamma_2}^\infty (1- \frac{1}{1+t \ x^{-1}})\ h(x)\ dx\right) $$
Now assume that $h(x)$ is a very complex expression,  I tell you that $\mathcal{L(t)}$ can not be integrated in closed form analytically and should be solved numerically.
So to solve my problem I have to take numerical values of $t$ in order to integrate numerically. But notice that to solve my main integration $F$ defined above, I need to set $$t=\sum_{i=1}^2\gamma_i^{-1} $$ and then integrate numerically to obtain $L\bigl(\sum_{i=1}^2\gamma_i^{-1} \bigl)$ by taking numerical values of $\gamma_1$ and $\gamma_2$.
Now that I took numerical values of $\gamma_1$ and $\gamma_2$ to solve my $\mathcal {L}(t)$. My inner integral is then numerical values for certain values of $\gamma_1$ and $\gamma_2$ that I picked. 
My question is, since I was obliged to take certain values of $\gamma_1$ and $\gamma_2$ to solve for my non-integrable function $\mathcal{L}(t)$ - since it happened to be that by $t$ is a function of $\gamma_1$ and $\gamma_2$. How would I proceed to solve my function $F$? Am I now obliged to integrate numerically over the same values of $\gamma_1$ and $\gamma_2$ I picked to solve for  $\mathcal{L}(t)$ ? 
In general do you think this way of analyzing such a problem is correct. Any better ideas on tackling such a problem? for example approximations/
Thanks
 A: Since the inner integral included in $\mathcal{L(t)}$ cannot be determined in closed form and must be solved by numerical methods, the double outer integral included in $F$ has to be solved numerically as well. Calculating the inner integral allows you to obtain a tabulation of the values of $\mathcal{L(t)}$ for several values of $\gamma_1$ and $\gamma_2$. From this tabulation, a 
numerical integration technique can be used to calculate $F$ as well.
Methods commonly used in this context (e.g. the Newton-Cotes quadrature formulas) typically approximate the integral of a function $f(x)$ by estimating various degree polynomials, starting from numerical values of $x$ and $f(x)$ tabulated at regularly spaced intervals. As correctly noted in the OP, it is therefore necessary to integrate numerically over the same values of $\gamma_1$ and $\gamma_2$ picked to solve for  $\mathcal{L}(t)$. In this view, it is convenient to pick $\gamma_1$ and $\gamma_2$ values by choosing regular intervals, whose dimension should be predefined according to the desired precision level.
To solve the problem of integrating on an infinite domain, you could consider the possibility of mapping the infinite interval to a finite interval using an appropriate substitution, i.e. choosing a function of $\gamma_1$ and $\gamma_2$ that goes to infinity in a finite range (for example $\tan\gamma_i$, $\cot \gamma_i$, or others).
