Basically you have to do it numerically!

But it is not very straightforward as you have singularities and \infty in limits. You have to resolve it or simply it analytically near the singularities (I say mathematically manipulate it in a way you can feed it to Wolfram Mathematica). Follow these steps:

1- Draw the numerator and denominator graphs for your range and parameters of interest to visually see singularities and behavior in infinity. (Remember you may need to break it down for your range of parameters if its behavior regime changes! like passing a bifurcation points based on tuning the parameters)

2- Try to use perturbation techniques to replace infinity related parts via symptomatically equivalent function (See the book by Van Dyke:

Perturbation Methods in Fluid Mechanics).

3- Either solve singular parts by analytically manipulating the integrant, or breaking its limit parts you can handle analytically and the parts you cannot. And finally check imaginary analysis techniques if possible to handle that. (See chapter 6 of the book: Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory. By Bender and Orszag)

4- Code parts 2 and 3 above, then numerically (with a mesh adoptive multiple integral solver) solve the integral while parts 2 and 3 are fed by a subroutines as a number. Based on part 1 you should know where singularities happen and you already went around them so you solver does not break.

5- If you want to verify your method is not wrong choose some dummy numbers and parameters which you have the analytical solution OR use you physical knowledge of the phenomenon. For example you know for this and this parameter you have to see this behavior from the physics

Hope this helps you!