Compute the (cumulative) probability function of $U = Y_1 - Y_2$, given that the joint density probability function $f(y_1, y_2)$ is
$$ \begin{align} f(y_1, y_2) &= e^{-y_1} \qquad \text{in} \qquad 0 \leq y_2 \leq y_1 \leq \infty \\ f(y_1, y_2) &= 0 \qquad \qquad \qquad \text{elsewhere} \end{align} $$
I get different results depending if I first integrate over $Y_1$, or over $Y_2$:
- $Y_1$-inner integration: $$ F_a(u) = P(Y_1 - Y_2 \leq u) = P(Y_1 \leq Y_2 + u) = \int_{0}^{\infty} \left(\int_{y_2}^{y_2+u} e^{- y_1} d y_1 \right) d y_2 = 1 - e^{-u} $$
- $Y_2$-inner integration: $$ F_b(u) = P(Y_2 \geq Y_1 - u) = \int_{u}^{\infty} \left( \int_{y_1 - u}^{y_1} e^{- y_1} d y_2 \right) d y_1 = u e^{-u} $$
How to correct the integration boundaries of one of the above, so that $F_a = F_b$?
$F_a(u)$ is correct. For $F_b(u)$, you missed another term, which is
$\int_0^u\left(\int_0^{y_1}e^{-y_1} dy_2\right) dy_1 = 1-e^{-u}-ue^{-u}$.