How to compute the integral:
$\int\limits_R (1+e^{-x_1}+e^{-x_2})^{-1}\ dx_2$
I tried substituting $e^{-x_2}=t$ simpliyfing thus to:-
$\int\limits_0^{\infty} (at+t^2)^{-1}\ dt$
where $a$ is a constant $=1+e^{-x_1}$.
The above can simply be integrated by the method of partial functions but it gives absurd results.
Answer is: $\frac{1}{a}$
I think $(1+e^{-x_1}+e^{-x_2})^{-1}$ may not be the joint density function but the joint cumulative distribution function $F(x_1,x_2) = \mathbb P(X_1 \le x_1, X_2 \le x_2)$, close to $0$ when either $x_1$ or $x_2$ are very negative and close to $1$ when they are both very positive
If you were to integrate it with respect to $x_2$, the indefinite integral would be $\frac{\log_e(e^{x_2}(1 + e^{-x_1})+1)}{1 + e^{-x_1}} + c$ which gives an unbounded result as $x_2$ increases
The joint density seems to be $\frac{2 e^{-x_1-x_2}}{(1 + e^{-x_1}+e^{-x_2})^3}$. You could integrate this over $x_2$ in $\mathbb R$ but there is an easier approach
To find the marginal cumulative distribution function $\mathbb P(X_1 \le x_1)$ you could simply let $x_2 \to +\infty$ to give $\mathbb P(X_1 \le x_1)=(1+e^{-x_1})^{-1}$ (your $\frac1a$) with a logistic marginal density $\frac{e^{-x_1}}{(1 + e^{-x_1})^2}$