ProblemCompute the flux of the vector field 
> F := x^2*i+x*y*j+x^3*y^3*k;
2 3 3 F := x i + x y j + x y k
over the surface of the tetrahedron bounded by the plane 

Recall the expression of the Divergence Theorem: 
> ;
/ / / / /        F.N dS =    div F dV      / / / / / S V
Let's start with the right hand side tripple integral. First we compute div F. 
> divF := diverge(evalm(F),[x,y,z]);
divF := 3 x
and we then integrate over the solid tetrahedron, 
> I2 := Int(Int(Int(divF,z=0..1xy),y=0..1x),x=0..1);
1 1  x 1  x  y / / /    I2 :=    3 x dz dy dx    / / / 0 0 0> I2 := int(int(int(divF,z=0..1xy),y=0..1x),x=0..1);
I2 := 1/8
Now the left hand side. Let us split the surface into the four triangles: S = Sxy + Sxz + Syz + Sxyz. Where Sab is the triangle on the ab plane and Sxyz is the triangle on the plane x+y+z=1. Let us denote the surface integrals over Sxy, etc... by J1,J2,J3 and J4. We have, that the outward unit normal for Sxy is (k) and thus, 
> ;
/ /   3 3 J1 :=   x y dS   / / Sxy
which evaluates to: 
> J1 := int(int(dotprod(F,k),y=0..1x),x=0..1);
1 J1 :=  1120
Now compute the two central surface integrals J2 and J3. The outward unit normals for Sxz and Syz are j and i respectively so, 
> ;
/ /   J2 :=   x y dS = 0   / / Sxz
this integral is 0 since y=0 on the xz plane. For J3 we have, 
> J3 := Int(Int(dotprod(F,i),S),u);
/ /   2 J3 :=   x dS = 0   / / Syz
again this integral is 0 since x=0 on the yz plane. For the last integral we see that the ouward unit normal is 
> N := evalm(vector([1,1,1])/sqrt(3));
[ 1/2 1/2 1/2] N := [1/3 3 , 1/3 3 , 1/3 3 ]
Notice that we need to divide by sqrt(3) (the norm of (1,1,1)) to get a UNIT normal. Hence, 
> ;
/ /   2 1/2 1/2 3 3 1/2 J4 :=   1/3 x 3 + 1/3 x y 3 + 1/3 x y 3 dS   / / Sxyz
Computing the dS for this last integral we find

> J4 := Int(Int(dotprod(F,[1,1,1]),y=0..1x),x=0..1);
1 1  x / /   2 3 3 J4 :=   x + x y + x y dy dx   / / 0 0
which evaluates to: 
> J4 := int(int(dotprod(F,[1,1,1]),y=0..1x),x=0..1);
141 J4 :=  1120
Finally combining the J's we get: 
> I1 := J1+0+0+J4;
I1 := 1/8
The same as before!! 