<s>[INST] This double integral $$  \int_0^1\int_0^1x^3y^2\sqrt{1+x^2+y^2}\hspace{1mm}dxdy$$
We have to compute this up to 4 decimal places
 [/INST] If you just need an approximation of this integral, you could use numerical methods. Matlab or Wolfram's Double Integral Calculator, for example, can do it just fine : 
$$ \int_{0}^{1} \int_{0}^{1} x^{3}y^{2}\sqrt{1+x^{2}+y^{2}} \, \mathrm{d}x \mathrm{d}y \simeq 0.125065 $$
But if you want an exact computation, this might be harder.
To begin with, I wouldn't use polar coordinates here because you integrate on the unit square $[0,1] \times [0,1]$ and the parametrization of this unit square in polar coordinates is not straightforward.
Let 
$$ I = \int_{0}^{1} \int_{0}^{1} x^{3}y^{2} \sqrt{1+x^{2}+y^{2}} \, \mathrm{d} x \mathrm{d}y.$$
Using Fubini's theorem, you can write :
$$I = \int_{0}^{1} y^{2} \Bigg( \int_{0}^{1} x^{3} \sqrt{1+y^{2}+x^{2}} \, \mathrm{d}x \Bigg) \, \mathrm{d}y. \tag{1}$$
Now, we have to compute an integral of this form :
$$ J=\int_{0}^{1} x^{3} \sqrt{a+x^{2}} \, \mathrm{d}x $$
with $a \in \mathbb{R}^{+}$. This can be done using two changes of variable :
$$
\begin{align*}
J &= {} \int_{0}^{1} x^{3} \sqrt{a+x^{2}} \, \mathrm{d}x \\[2mm]
 & \mathop{=} \limits_{u=x^2} \frac{1}{2} \int_{0}^{1} u \sqrt{a+u} \, \mathrm{d}u \\[2mm]
 & = \frac{1}{2} \int_{a}^{a+1} (v-a)\sqrt{v} \, \mathrm{d}v \\[2mm]
 & = \frac{1}{5} (a+1)^{\frac{5}{2}} + \frac{2}{15}a^{\frac{5}{2}}-\frac{a}{3}(a+1)^{\frac{3}{2}} \\
\end{align*}
$$
As a consequence, $(1)$ becomes :
$$ I = \int_{0}^{1} y^{2} \Big[ \frac{1}{5}(2+y^{2})^{\frac{5}{2}} + \frac{2}{15}(1+y^{2})^{\frac{5}{2}} - \frac{1}{3}(1+y^{2})(2+y^{2})^{\frac{3}{2}} \Big] \, \mathrm{d}y. $$
At this point, we have to compute the three following integrals : 
$$I_{1} = \int y^{2}(1+y^{2})^{\frac{5}{2}} \, \mathrm{d}y, \quad I_{2} = \int y^{2}(1+y^{2})^{\frac{3}{2}} \, \mathrm{d}y \quad \mathrm{and} \quad I_{3} =\int (1+y^{2})^{\frac{3}{2}} \, \mathrm{d}y.$$
One method to compute these integrals (it might not be the easiest way) would be to consider the change of variable $y = \sinh(t)$. Even though, the computation is not completely straightforward. If you want, I can give more explanations about the computation of these integrals. Wolfram gives the following results :
$$
\begin{align*}
I_{1} & = {} \frac{1}{8} \Big( 19\sqrt{3} - 5 \sinh^{-1}\big(\frac{1}{\sqrt{2}} \big) \Big) \\
I_{2} & = \frac{1}{384} \big( 317\sqrt{2} - 15 \sinh^{-1}(1) \big) \\
I_{3} & = \frac{1}{8} \Big( 11\sqrt{3} - \sinh^{-1}\big( \frac{1}{\sqrt{2}} \big) \Big)
\end{align*}
$$
In the end :

$$I = \frac{317\sqrt{2} - 15 \sinh^{-1}(1)}{2880} + \frac{1}{40}\Bigg(19\sqrt{3}-5\sinh^{-1} \Big( \frac{1}{\sqrt{2}} \Big) \Bigg)+\frac{1}{24}\Bigg( \sinh^{-1}\Bigg( \frac{1}{\sqrt{2}} \Big) - 11\sqrt{3} \Bigg).$$

</s>