Given the joint density $f(x,y)$ we are asked to find $P(X + Y \leq c)$ where $c$ is some constant non-negative number.
My reasoning for this problem is thus:
$P(X+Y \leq c) = P(X \leq c - Y)$
In order to find this probability, we need the marginal probability distribution: $P(X)$. We can then integrate $P(x)$ from the lower bound (say, $0$) to $c-y$ for the desired probability.
$P(X)$ is found by taking the integral of the joint density with respect to y from the lower bound to the upper bound (say $0$ to $1$).
$P(X) = \int_{0}^{1} f(x,y)dy$
Then we just take the integral as described above and we have the following double integral:
(1) $$\int_{0}^{c - y} \int_{0}^{1} f(x,y)dydx$$
However, in all the examples I've seen thus far, the integrals are reversed:
(2) $$\int_{0}^{1} \int_{0}^{c-y} f(x,y)dxdy$$
Doing some research I discovered Fubini's theorem which at first seemed to explain why this switch is possible (if not explain why it's necessary).
However, I soon discovered that the two above integrals are not in fact equal (e.g if $f(x,y) = \frac {6}{7} (x+y)^{2}$ ). Which means that my reasoning is wrong, Fubini's theorem doesn't apply here, and I should arrive at the second integral from the very beginning.
Can some explain where my reasoning is faulty and how one can arrive at the second integral?
Generally, for continuous real-valued random variables, $X$ and $Y$ with a joint probability density function $f_{X,Y}$, Fubini's theorem allows that: $$\begin{align}\mathsf P(X+Y\leqslant c) ~&= \int_{\{(x,y):x+y\leqslant c\}} f_{X,Y}(x,y)~\mathrm d(x,y) \\[1ex] &= \int_{-\infty}^{\infty}\int_{-\infty}^{c-x} f_{X,Y}(x,y)~\mathrm d y~\mathrm d x \\[1ex] &= \int_{-\infty}^{\infty}\int_{-\infty}^{c-y} f_{X,Y}(x,y)~\mathrm d x~\mathrm d y\end{align}$$
(Notice that the constraint is placed on the inner integral and is a function of the bound variable for the outer integral.)
In the specific case when the supports for the probability density functions of $X,Y$ are strictly non-negative, and $c$ is too, this may be expressed as:
$$\begin{align}\mathsf P(X+Y\leqslant c) ~&= \int_{\{(x,y):0\leqslant x+y\leqslant c\}} f_{X,Y}(x,y)~\mathrm d(x,y) \\[1ex] &= \int_{0}^{c}\int_{0}^{c-x} f_{X,Y}(x,y)~\mathrm d y~\mathrm d x \\[1ex] &= \int_{0}^{c}\int_{0}^{c-y} f_{X,Y}(x,y)~\mathrm d x~\mathrm d y\end{align}$$
(Notice the constraint on the outer integral is not a function of either bound variables.)
In the particular case when $f_{X,Y}(x,y) = \tfrac 67(x+y)^2~\mathbf 1_{0\leqslant x\leqslant 1, 0\leqslant y\leqslant 1}$, we have:
$$\begin{align}\mathsf P(X+Y\leqslant c) ~&=~ \int_{\{(x,y):0\leqslant x\leqslant 1,0\leqslant y\leqslant 1,0\leqslant x+y\leqslant\min\{2,c\}\}}\tfrac 67(x+y)^2~\mathrm d (x,y) \\[2ex] &=~ \int_0^{\min\{1,c\}}\int_0^{\min\{1,c-x\}} \tfrac 67 (x+y)^2~\mathrm d y~\mathrm d x \\[2ex]&=~ \begin{cases}0&:& c<0 \\ \int_{0}^{c}\int_{0}^{c-x} \tfrac 67(x+y)^2~\mathrm d y~\mathrm d x &:& 0\leqslant c < 1\\ \int_{0}^{c-1}\int_{0}^{1} \tfrac 67(x+y)^2~\mathrm d y~\mathrm d x+\int_{c-1}^{1}\int_{0}^{c-x} \tfrac 67(x+y)^2~\mathrm d y~\mathrm d x &:& 1\leqslant c < 2\\ 1 &:& 2\leqslant c \end{cases}\\[2ex] &= \begin{cases}0&:& c<0 \\ \int_{0}^{c}\int_{0}^{c-y} \tfrac 67(x+y)^2~\mathrm d x~\mathrm d y &:& 0\leqslant c < 1\\ \int_{0}^{c-1}\int_{0}^{1} \tfrac 67(x+y)^2~\mathrm d x~\mathrm d y + \int_{c-1}^{1}\int_{0}^{c-y} \tfrac 67(x+y)^2~\mathrm d x~\mathrm d y&:& 1\leqslant c < 2\\ 1 &:& 2\leqslant c \end{cases}\end{align}$$
Evaluation is left to you.