The functions are given on the $xy$ plane, the $xz$ plane and the $yz$ plane, respectively, i.e. $y=1-x^2$, $z=1-x^2$ and $z=1-y^2$ for $0 \leq x \leq 1, \ 0 \leq y \leq 1$ and $0 \leq z \leq 1$.
My original purpose was to find the area bounded by the three functions and the lines $x=0, y=0$ and $z=0$.
To do so I tried to find the function $z=f(x, y)$ that allows me to calculate
$$V=\int_{x=0}^{x=1} \int_{y=0}^{y=1-x^2} f(x, y) \ dy \ dx$$
but I couldn't. The best I could find was $f(x,y)=1-x^2-y^2$ which obviously does not match with $y=1-x^2$ when $z=0$.
How can I find it?
The region $\mathcal{R}$ defined by the given constraints is the intersection of three solids. The boundary surfaces of each solid are a parabolic cylinder and four planes. You can only solve this problem by using two different functions $z=f(x,y)$ and $z=g(x,y)$. There are only two possibilities: $z=1-x^2$ and $z=1-y^2$.
The intersection of the surfaces $z=1-x^2$ and $z=1-y^2$ is the plane $y=x$. The plane $y=x$ intersects the surface $y=1-x^2$ at $x=(\sqrt{5}-1)/2 = \phi$. Define the points $A,B,C,D$ in the $xy$ plane as $A=(0,0)$, $B=(1,0)$, $C=(0,1)$, and $D=(\phi,\phi)$. The $xy$ domain of the region $\mathcal{R}$ when $z=1-x^2$ is the triangle $T_0$ with the two linear edges $AB$ and $AD$ and the curved edge $y=1-x^2$ for $x\in [\phi,1]$. We could integrate over $T_0$ with one double integral if we integrate with respect to $x$ first, but this would result in some square roots in the integrals. So we break $T_0$ into two triangles $T_1$ and $T_2$. The first triangle $T_1$ is the right triangle with vertices $A$, $D$, and $E=(\phi,0)$. The second triangle $T_2$ is the triangle with linear edges $EB$ and $ED$ and the curved edge of $T_0$. Computing the double integrals over $T_1$ and $T_2$ is easy: $$ \begin{align} I_1 &= \iint_{T_1} 1-x^2 \,dA = \int_0^{\phi}(1-x^2)\int_0^x \,dy\,dx \\ &=\int_0^{\phi} x(1-x^2) \, dx = \frac{\phi^2}{2} - \frac{\phi^4}{4}. \end{align} $$ and $$ \begin{align} I_2 &= \iint_{T_2} 1-x^2 \,dA = \int_{\phi}^1 (1-x^2) \int_0^{1-x^2} dy\,dx \\ &= \int_{\phi}^1 (1-x^2)^2 \,dx = \frac{8}{15} - \phi + \frac{2}{3}\phi^3 - \frac{1}{5}\phi^5. \end{align} $$
The $xy$ domain of the region $\mathcal{R}$ when $z=1-y^2$ is the triangle $T_3$ with linear edges $AC$ and $AD$ and curved edge $y=1-x^2$ for $x\in [0,\phi]$. The double integral over $T_3$ is $$ \begin{align} I_3 &= \iint_{T_3} 1-y^2 \,dA = \int_0^{\phi}\int_x^{1-x^2} 1-y^2 \,dy\,dx \\ &= \frac{2}{3}\phi - \frac{1}{2}\phi^2 + \frac{1}{12}\phi^4 - \frac{1}{5}\phi^5 + \frac{1}{21}\phi^7. \end{align} $$ The total volume is $$ \begin{align} V &= I_1+I_2+I_3 = \frac{8}{15} - \frac{1}{3}\phi + \frac{2}{3}\phi^3 - \frac{1}{6}\phi^4 - \frac{2}{5}\phi^5 + \frac{1}{21}\phi^7 \\ &= \frac{123}{420} + \frac{5}{84}\sqrt{5} \approx 0.425956. \end{align} $$
One can easily check this result using Monte Carlo integration with $N$ uniform random samples selected from the unit cube. For each sample inside the region $\mathcal{R}$, add 1 to a variable $S$ that is initialized to 0. The Monte Carlo estimate is $S/N$.