Systematic method to change the order of integration in multiple integrals In many examples of computation of multiple integrals, it is necessary to change the order of integration to achieve the computation. For example, $I=\int_0^1\int_y^1 \cos(x^2)\ dx\ dy$ can be computed using Fubini's theorem as $I=\int_0^1\int_0^x\cos(x^2)\ dy\ dx=\int_0^1 x\cos(x^2)\ dx=\left[\frac{cos(x^2)}{2}\right]_0^1=\frac{\cos(1)}{2}$. 
It is generally easy to change the order of integration in double integral, but way more difficult in some cases of triple (or more) integrals. 
For example, 
$\begin{array}{rcl}\int_0^1\int_0^{1-x^2}\int_0^{1-x}f(x,y,z)\ dy\ dz\ dx &
= & 
\int_0^1\int_0^{1-\sqrt{1-z}}\int_0^{\sqrt{1-z}}f(x,y,z)\ dx\ dy\ dz  \\ & + &  \int_0^1\int_{1-\sqrt{1-z}}^{1}\int_0^{1-y}f(x,y,z)\ dx\ dy\ dz
\end{array}$ 
Is there a general process to obtain this result? Please note that I know how to tackle this example, I am looking for a generic method/algorithm (hopefully simple enough to explain it to high school students.) Thank you.
 A: In the first place, you should disambiguate the integration bounds, indicating to what variable they relate, and show the nesting, like
$\int_{y=0}^{1}[\int_{x=y}^{1}\cos(x^2)\ dx]\ dy=\int_{x=0}^{1}[\int_{y=0}^{x}\cos(x^2)\ dy]\ dx$.
Now, the domain of integration is represented by the inequations
$0\le y\le1$ and $y\le x\le 1$, a triangle inscribed in the tile $[0,1]\times[0,1]$. 
Choosing a value for $y$, let $Y$, gives you a range for $x$: $Y\le x\le 1$. Similarly, choosing a value for $x$ gives you a range for $y$: $0\le y\le X$.
Your second example is presumably based on the inequations $0\le x\le 1$, $0\le z \le 1-x$ and $0\le y\le 1-x^2$, a truncated parabolic cone inside $[0,1]\times[0,1]\times[0,1]$.
Choosing a value for $z$ gives you a domain for $(x, y)$: $0\le Z\le 1-x$, i.e. $0\le x\le Z$ or $1-Z\le x\le 1$, and $0\le y \le 1-x^2$.
Then choosing a value for $y$ gives you a range for $x$: $x^2\le 1-y$, i.e. $-\sqrt{1-y}\le x\le \sqrt{1-y}$.
Hence: $\int_{z=0}^1[\int_{y=0}^{z}[\int_{x=-\sqrt{1-z}}^{\sqrt{1-z}}dx]dy+\int_{y=1-z}^{1}[\int_{x=-\sqrt{1-z}}^{\sqrt{1-z}}dx]dy]dz$
The general methodology is to write down the inequations describing the integration domain, then to choose the first integration variable and to "slice" the domain by making this variable a constant. You will obtain a sub-domain in $d-1$ dimensions, and you can repeat the process by choosing the next independent integration variable.
