In this case you don't need Fubini's theorem or equivalent to produce an iterated integral (although a careless exposition can conceal this fact), since what one is actually doing is using the linearity of the integral twice: you have
$$ J^2 = J\int_0^{\infty} e^{-x^2} \, dx = \int_0^{\infty} e^{-x^2} J \, dx. $$
Now, $e^{-x^2}$ is also a constant when integrating with respect to $y$, so we have
$$ e^{-x^2} J = e^{-x^2} \int_0^{\infty} e^{-y^2}\, dy = \int_0^{\infty} e^{-x^2}e^{-y^2} \, dy, $$
and so
$$ J^2 = \int_0^{\infty} \left( \int_0^{\infty} e^{-(x^2+y^2)} dy \right) dx $$
as desired.
The inner integral $\int_0^{\infty} e^{-(x^2+y^2)} dy$ is now a function of $x$, and we substitute to write this function in a different way, as $x\int_0^{\infty} e^{-x^2(1+t^2)} \, dt $. Now is when you are forced to change the order of integration using some theorem or other. Either we embrace a Lebesgue integral version, such as Tonelli's theorem,
Let $f : (a,b) \times (c,d) \to \mathbb{R}$ be positive and Lebesgue-integrable, where $-\infty \leq a,b \leq \infty$ and $ -\infty \leq c,d \leq \infty$. Then
$$ \int_{(a,b)\times (c,d)} f(x,y) \, dx \times dy = \int_a^b \left( \int_c^d f(x,y) \, dy \right) dx = \int_c^d \left( \int_a^b f(x,y) \, dx \right) dy. $$
Here, though, we can manage with a simple version for Riemann integrals:
Let $f : (a,b) \times (c,d) \to \mathbb{R}$ be positive and Riemann-integrable, where $-\infty <a,b < \infty$ and $ -\infty < c,d < \infty$. Then
$$ \int_{(a,b)\times (c,d)} f(x,y) \, dx \times dy = \int_a^b \left( \int_c^d f(x,y) \, dy \right) dx = \int_c^d \left( \int_a^b f(x,y) \, dx \right) dy. $$
In particular, $x = 1-(1-u)^{-1}$ provides a bijection between $(0,1)$ and $(0,\infty)$, so we find that
$$ J^2 = \int_0^1 \int_0^1 \frac{u}{(1-u)^3(1-v)^2}\exp{\left( - \frac{u^2}{(1-u)^2}\left( 1 + \frac{v^2}{(1-v)^2} \right) \right)} \, dv \, du. $$
It turns out that the integrand is in fact continuous on the whole open square, so it's integrable, and we can use the lightweight Tonelli to swap the integrals, and then after undoing the substitutions, we're finally at
$$ J^2 = \int_0^{\infty} \int_0^{\infty} xe^{-x^2(1+t^2)} \, dx \, dt, $$
which is straightforward to do. Alternatively if one is happy to use the Lebesgue integral, we don't need to muck about with substitutions, and can apply the proper version of Tonelli to the integral on $(0,\infty) \times (0,\infty)$.
Best Answer
$$ \int _0^1\int _0^xf(x,y)dydx=\int _0^1\int _0^1\chi _{[0,x]}(y)f(x,y)dydx $$
Now apply Fubini to get
$$ \int _0^1\int _0^xf(x,y)dydx=\int _0^1\int _0^1\chi _{[0,x]}(y)f(x,y)dxdy=\int _0^1\int _y^1f(x,y)dxdy, $$
where I have used the fact that $\chi _{[0,x]}(y)=0$ unless $y\leq x\leq 1$.
Techincally speaking, you can only apply Fubini (or Tonelli) for a rectangular region. To do more general regions, you have to play around with characteristic functions as I just did (or possibly even do a change of variables) and then apply Fubini (or Tonelli). However, in practice, with a bit of geometric intuition, you can figure out what the bounds should be without doing this.