How to handle complicated double integral like this numerically

Question

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?

Answer 1

Thanks to the help from David, I managed to solve it using integral2 function in Matlab:

c=0.13;
b1=8.64938;
b2=11.1234
q1=22.07441;
q2=32.3352;
fun=@(t,y)(y-t).*b1.*b2.*c.^2.*t.^(c-1).*y.^(c-1).*exp(q1+q2-b1.*t.^c-b2.*y.^c)./((1+exp(q1-b1.*t.^c).^2.*(1+exp(q2-b2.*y.^c).^2)))
syms t y
tmax=@(y)y
integral2(fun,0,+inf,0,tmax)