This pdf seems related to the Weibull distribution. I made some abrreviations
Etta_d = xid Etta_D = xiD
Gamma_d = gd Gamma_D = gD
Mu_d = mud Mu_D = muD,
and changed the sign of the argument of the second exp, multiplying by it rather than dividing. This leads to
PDF = @(x,y) (gd*gD)/(Del*xid*xiD)*(1./y) ...
.*exp(-(1 -mud*gd +gd*x./y ).^(-1/xid)).*(1 -mud*gd +gd*x./y ).^(-1-1/xid) ...
.*exp(-(1 -muD*gD +gD*y/Del).^(-1/xiD)).*(1 -muD*gD +gD*y/Del).^(-1-1/xiD); 
Take a look at line  , considered as a distribution on its own. There is a basic normalized distribution (whose name I don't know),
(1/xiD) * exp(-y.^(-1/xiD)) .* y.^(-1-1/xiD)
with 0 <= y <= inf. Variable y cannot be negative since it is taken to a fractional power.
Line  contains a scaled and shifted version of y. Now the argument itself must be positive,
0 <= (1 -muD*gD +gD*y/Del) <= inf
The net result is that y can take on some negative values, and the mean value of y turns out to be y_mean = 1, which I assume is intentional. The lower limit for y is
which is negative. If you integrate line  from y0 to inf, and include the normalization factors in front that are associated with that line,
you get 1.
For the full pdf, integration over y is complcated because now you have the condition in line 
0 <= (1 -mud*gd +gd*x./y ) <= inf
and the y limits depend on x. But you sensibly wanted to do both the x integration and the y integration as a check, which is fairly easy. Substitute
x = y*z dx = y*dz
and do the z integration first. The y in y*dz cancels the 1/y in the pdf and the resulting integration is from z0 to inf, where z0 can be calculated in the same way that y0 was. The result is
Del = 25;
xid = .02;
xiD = .08;
mud = -2.95;
muD = .05;
gD = (gamma(1-xiD)-1)/(1-muD);
gd = (gamma(1-xid)-1)/(1-mud);
PDF1 = @(z,y) (gd*gD)/(Del*xid*xiD) ...
.*exp(-(1 -mud*gd +gd*z ).^(-1/xid)).*(1 -mud*gd +gd*z ).^(-1-1/xid) ...
.*exp(-(1 -muD*gD +gD*y/Del).^(-1/xiD)).*(1 -muD*gD +gD*y/Del).^(-1-1/xiD);
z0 = ((mud*gd-1)/gd)*.9999;
y0 = ((muD*gD-1)*Del/gD)*.9999;
CDF1 = integral2(PDF1,z0,inf,y0,inf)
There are some integration issues near z0 and y0. To address that, I just multiplied each limit by .9999. To really do it right would take some work but I think the simple approach here shows that the result is correct.