Finding a function $z=f(x, y)$ for the double integral 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$.

