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
The region of integration is given by the inequalities $0\le x\le1$, $0\le y\le1-x$ and $0\le z\le1-x^2$. The latter two are equivalent to $x\le 1-y\le 1$ and $x^2\le1-z\le1$ (which is equivalent to $x\le\sqrt{1-z}\le1$).
Your integral is over the region defined by $0\le y\le1$, $0\le z\le1$ and $0\le x \le\sqrt{1-z}$. This is not the same as the original region. There are points like $(x,y,z)=(3/5,1,16/25)$ in your region which are not in the original region. You have neglected the condition that $0\le y\le 1-x$.
As you say, we must have $0\le y\le 1$ and $0\le z\le1$. Then the condition on $x$ is that $0\le x\le\min(1-y,\sqrt{1-z})$. When $0\le y\le1-\sqrt{1-z}$ that amounts to $0\le x\le \sqrt{1-z}$ and when $1-\sqrt{1-z}\le y\le 1$ it amounts to to $0\le x\le1-y$. So we get Stewart's limits.
By the way, is Stewart's book really 1,950+ pages?