By the divergence theorem, the surface integral given by
[tex]\displaystyle\iint_{\partial S}\mathbf F\cdot\mathbf n\,\mathrm dS[/tex]
(where the integral is computed over the entire boundary of the surface) is equivalent to the triple integral
[tex]\displaystyle\iiint_R\nabla\cdot\mathbf F\,\mathrm dV[/tex]
where [tex]V[/tex] is the volume of the region [tex]R[/tex] bounded by [tex]\partial S[/tex].
You have
[tex]\mathbf F(x,y,z)=x^2y\,\mathbf i+xy^2\,\mathbf j+4xyz\,\mathbf k[/tex]
[tex]\implies \nabla\cdot\mathbf F=\dfrac{\partial}{\partial x}[x^2y]+\dfrac{\partial}{\partial y}[xy^2]+\dfrac{\partial}{\partial z}[4xyz]=8xy[/tex]
and so the integral reduces to
[tex]\displaystyle8\int_{x=0}^{x=5}\int_{y=0}^{y=1}\int_{z=0}^{z=5-x-5y}xy\,\mathrm dz\,\mathrm dy\,\mathrm dx=-\dfrac{250}3[/tex]
