[tex]\vec F(x,y,z)=(x^2,-y^2,z^2)[/tex] has divergence
[tex]\mathrm{div}\vec F(x,y,z)=2x-2y+2z[/tex]
so that by the divergence theorem, the flux of [tex]\vec F[/tex] across the boundary of [tex]D[/tex] is
[tex]\displaystyle\iint_{\partial D}\vec F\cdot\mathrm d\vec S=2\iiint_D(x-y+z)\,\mathrm dV[/tex]
which you can compute by subtracting [the integral of [tex]\mathrm{div}\vec F[/tex] over the region bounded by [tex]9-x-y[/tex] and the coordinate axes] and [the integral of [tex]\mathrm{div}\vec F[/tex] over the region bounded by [tex]3-x-y[/tex] and the coordinate axes].
[tex]=2\displaystyle\left\{\int_0^9\int_0^{9-x}\int_0^{9-x-y}-\int_0^3\int_0^{3-x}\int_0^{3-x-y}\right\}(x-y+z)\,\mathrm dz\,\mathrm dy\,\mathrm dx[/tex]
[tex]=2\left(\dfrac{2187}8-\dfrac{27}8\right)=\boxed{540}[/tex]
