I have an improper integral which does not have a closed-form expression.
$$p1 = 2\pi \lambda_M \int_{r=0}^{\infty} \frac{e^{- \mathcal{Z} r^2}}{1 + \mathcal{Y}r^\alpha} r \, dr$$
If I used numerical-integration using simpson's method, how can I find the ideal value of the upper limit of $r$. I am having different results for different values of $r$
On
Sub $r=\sqrt{x/Z}$ and the integral becomes
$$\frac{\pi \lambda_M}{\sqrt{Z}} \int_0^{\infty} dx \, \frac{e^{-x}}{1+Y Z^{-\alpha/2} x^{\alpha/2}} $$
Now you can use something called Gauss-Laguerre quadrature. This is where one uses the orthogonality properties of Laguerre polynomials $L_n$ to form the approximation
$$\int_0^{\infty} dx \, f(x) e^{-x} \approx \sum_{i=1}^n w_i f(x_i)$$
where
$$L_n(x_i) = 0$$
and
$$w_i = \frac{x_i}{(n+1)^2 L_{n+1}(x_i)^2} $$
Note that the Laguerre polynomial is defined as
$$L_n(x) = e^x \frac{d^n}{dx^n} (e^{-x} x^n) $$
A table of zeros may be found here.
One way is to not use an upper truncation limit at all and instead map to an integral over $(0,1)$ using a substitution such as $x=1/(1+r)$, or $r=\tan(\pi x/2)$.
Another way is to note that the integrand is bounded by $e^{-\mathcal{Z} r^2}$ (assuming $\mathcal{Y} \geq 0$ so that everything makes sense). This combined with some familiar estimates for tails of the normal distribution allows you to show that the tail of your integral (without the coefficient in front) is bounded by $e^{-\mathcal{Z} R^2}$.
Still another way would be to use Gauss-Hermite quadrature (which would converge quite rapidly if $\mathcal{Y}$ isn't too large).