More generally, consider the distribution whose density is
$$f_p(u, \theta) = \theta C_p\int_1^\infty v^{-p} e^{-\theta u v}\,\mathrm{d}v = \theta C_pE_p(\theta u)$$
for $u \gt 0$ and parameters $p\gt 0,$ $\theta \gt 0.$ $E_p$ is the Exponential Integral function and $C_p$ is a normalizing constant (which we will find at the end). Notice, for future reference, that
$$E_{p+1}(0) = \int_1^\infty v^{-(p+1)}\,\mathrm{d}v = \frac{1}{p}.$$
Moreover, since for any $z\gt 0$
$$
0\le E_p(z) = \int_1^\infty v^{-p} e^{-zv}\,\mathrm{d}v\le \int_1^\infty e^{-zv}\,\mathrm{d}v = \frac{1}{z},$$
it follows from the Squeeze Theorem (for limits) that
$$\lim_{z\to\infty} E_p(z) = 0.$$
The distribution function (cdf) is, by definition,
$$F_p(u,\theta) = \int_0^u f_p(z)\,\mathrm{d}z = C_p \int_0^u E_p(\theta z)\,\theta\mathrm{d}z = C_p \int_0^{\theta u} E_p(z)\,\mathrm{d}z.$$
The integrand defining $E_p$ is smooth and decreasing so quickly at its upper limit that we may interchange the operations of differentiating and integrating, showing
$$\begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}z}E_{p+1}(z) &= \frac{\mathrm{d}}{\mathrm{d}z}\int_1^\infty v^{-(p+1)} e^{-zv}\,\mathrm{d}v\\&=\int_1^\infty \frac{\mathrm{d}}{\mathrm{d}z}v^{-(p+1)} e^{-zv}\,\mathrm{d}v\\&=-\int_1^\infty v^{-p} e^{-zv}\,\mathrm{d}v\\&=-E_p(z).
\end{aligned}$$
It is immediate (by the Fundamental Theorem of Calculus) that
$$F_p(u,\theta) = C_p \left[E_{p+1}(0) - E_{p+1}(\theta u)\right] = C_p \left[\frac{1}{p} - E_{p+1}(\theta u)\right].$$
Finally, since $\lim_{u\to\infty}F_p(u,\theta)=1$ (by the Law of Total Probability), we find
$$1 = \lim_{u\to\infty}C_p \left[\frac{1}{p} - E_{p+1}(\theta u)\right] = \frac{C_p}{p},$$
showing $C_p = p$ and giving
$$F_p(u,\theta) = 1 - pE_{p+1}(\theta u).$$
Set $p=5$ and $\theta=4$ for the answer to the question. Here are graphs of the density and its cdf:
These were drawn in gray using the analytical solution. As a check, over them are plotted, in red, the values obtained from numerically computing the single integral for $f$ and the double integral for $F.$ The R code that generated these graphs follows.
#
# Solution.
#
library(expint)
f <- function(u, p, theta) theta * p * expint_En(theta * u, p)
F <- function(u, p, theta) 1 - p * expint_En(theta * u, p+1)
#
# Brute force verification using numerical integration.
#
f. <- Vectorize(function(u, p, theta, ...) {
theta * p * integrate(function(v) exp(-theta * u * v) / v^p, 1, Inf, ...)$value
}, "u")
F. <- Vectorize(function(u, p, theta, ...) {
integrate(function(z) f.(z, p, theta), 0, u, ...)$value
}, "u")
#
# Example.
#
p <- 5
theta <- 4
par(mfrow=c(1,2))
curve(f(x, p, theta), 0, 1.5, n=501, lwd=2, col=gray(0.25),
xlab="u", ylab="Density",
main=bquote(f[.(p)](u, .(theta))))
curve(f.(x, p, theta), add = TRUE, col="Red", lwd=2, lty=2)
curve(F(x, p, theta), 0, 1.5, n=501, lwd=2, col=gray(0.25),
yaxp=c(0,1,1),
xlab="u", ylab="Probability",
main=bquote(F[.(p)](u, .(theta))))
curve(F.(x, p, theta), add = TRUE, col="Red", lwd=2, lty=2)
par(mfrow=c(1,1))
Best Answer
As indicated in the earlier comments, once you get a sample from the joint distribution of $(X_1,X_2,X_3)$, $$(x_1^1,x_2^1,x_3^1),\ldots,(x_1^t,x_2^t,x_3^t)$$ the marginal sample $$(x_1^1,x_2^1),\ldots,(x_1^t,x_2^t)$$is indeed a sample from the marginal joint distribution of $(X_1,X_2)$ and you can ignore the simulated $x_3^j$'s. They can however be useful in Monte Carlo evaluations through a technique called Rao-Blackwellisation since the average$$\frac{1}{t}\sum_{i=1}^t h(x_1^i,x_2^i)$$is improved by the average$$\frac{1}{t}\sum_{i=1}^t \mathbb{E}[h(x_1^i,x_2^i)|x_3^i]$$as a conditional expectation shares the same expectation as the original but reduces the variance. See for instance this discussion on Cross Validated.