# How to handle complicated double integral like this numerically

I'm now frustrated by this complicated double integral when doing my thesis work recently:

$$\int _0^{+\infty }\int _0^y\:\left(y-x\right)·\frac{b_2·c·y^{c-1}}{\left(1+e^{q_2{-b_2y^c}}\right)^2}·\frac{b_1·c·x^{c-1}}{\left(1+e^{q_1-b_1x^c}\right)^2}dxdy$$

where $b_1=8.64938, b_2=11.1234, q_1=22.07441, q_2=32.3352, c=0.13$.

I tried to solve it in Matlab but no go:

c=0.13;
b1=8.64938;
b2=11.1234
q1=22.07441;
q2=32.3352;
function=(y-x)*b1*b2*c^2*x^(c-1)*y^(c-1)*exp(q1+q2-b1*x^c-b2*y^c)/((1+exp(q1-b1*x^c)^2*(1+exp(q2-b2*y^c)^2)
syms x y
double(int(int(function,x,0,y),y,0,+inf))


I also looked into R and some packages like cubature, however, I couldn't find a package fitting my problem. So, is there a proper way to handle it numerically, by using R or Matlab or anything else?

• Look at Matlab's numerical integration functions integral2 (best choice), quadgk, trapz. – David Oct 11 '15 at 22:53

c=0.13;