This is example 8.9(a) in Rudin's Real and Complex Analysis, (alternatively, exercise 10.2 in Rudin's Principles of Mathematical Analysis).
Let $X$ and $Y$ be the closed unit interval $[0,1]$, let $\{\delta_n\}$ be an increasing sequence of distinct points in $[0,1]$ that converges to $1$, and to each positive integer $n$, let $g_n$ be a real continuous function on $[0,1]$ with support in $(\delta_n,\delta_{n+1})$, and such that $\int_{0}^{1} g_n(t)~dt=1$. Define $f$ over $X\times Y$ as follows: $$ f(x,y):=\sum_{n=1}^\infty[g_n(x)-g_{n+1}(x)]g_n(y). $$ If $x$ is held fixed, then $g_n(x)-g_{n+1}(x)$ is nonzero for at most two choices of $n$, and $f$ reduces to a finite sum. So, \begin{align*} \int_X~dx\int_Y f(x,y)~dy&=\int_0^1~dx\int_0^1\sum_{n=1}^\infty[g_n(x)-g_{n+1}(x)]g_n(y)~dy \\ &=\int_0^1~dx\sum_{n=1}^\infty\int_0^1[g_n(x)-g_{n+1}(x)]g_n(y)~dy \\ &=\int_0^1\sum_{n=1}^\infty g_n(x)-g_{n+1}(x)~dx \\ &=\int_0^1 g_1(x)~dx=1. \end{align*} On the other hand, if $y$ is held fixed, then $g_n(y)$ is nonzero for at most one choice of $n$, and $f$ reduces to a finite sum. So, \begin{align*} \int_Y~dy\int_X f(x,y)~dx&=\int_0^1~dy\int_0^1\sum_{n=1}^\infty[g_n(x)-g_{n+1}(x)]g_n(y)~dx \\ &=\int_0^1~dy\sum_{n=1}^\infty\int_0^1[g_n(x)-g_{n+1}(x)]g_n(y)~dx\\ &=\int_0^1 \sum_{n=1}^\infty g_n(y)-g_n(y)~dy\\ &=\int_0^1 0~dy = 0. \end{align*}
The double integrals are not interchangeable in this example, in part, because $$ \int_0^1~dx\int_0^1|f(x,y)|~dy=\infty. $$ Why is this true?