Evaluate the integral of the vector field $$F(x,y,z)=(xy,y^2,y^3)$$ over the boundary of the unit cube.
My method is to parametrize all the six boundaries of the unit cube in the following way $$P_1=(u,v,0),\quad P_2=(u,v,1), \quad P_3=(0,u,v)..., \quad 0\le u,v\le 1$$ $$P_4=(1,u,v), \quad P_5=(u,0,v), \quad P_6=(u,1,v).$$
then evaluate $F(P_i)$ which i will denote it as $F_i$, and then evaluate $$\frac{\partial P_i}{\partial u}\times\frac{\partial P_i}{\partial v}=N_i$$
and finally $$\iint_\limits{S}F\cdot n\,d\sigma=\sum_{i=1}^6\left(\iint_\limits{R}F_i\cdot N_i\,du\,dv\right)$$ It turns out that the region $R$ is just $[0,1]\times[0,1]$, now the problem is that I got $0$ as an answer but the book say it is $3/2$, so is it just a calculation error or my method doesn't work?
You can't parameterize the top and bottom surface as $(u,v,0)$ and $(u,v,1)$ -- that means that your normal vectors $N_1$ and $N_2$ will be equal, and this vector points into the cube at one surface and out of the cube at the other.
What you're doing ought to work once you fix the inconsistent directions of your normals. However, if you have the divergence theorem available, you're almost certainly expected to use that instead of integrating over each of the six sides separately. That will be much quicker, especially because the divergence of this particular field is much simpler than the field itself.