We can do this 1 of 2 ways. The first of which is to compute the flux integrals for the 3 distinct sides of the cylinder (top, bottem and side) or use the divergence theorem.
I will use the divergence theorem,
There should be a closed symbol on the right side but I dont know the latex code for this :P
Our Function, has the divergence such that,
Our flux is therefore,
Noting our dz bounds are the bounds we defined for z at the start of this problem, and our x and y bounds can be most easily expressed in cylindrical co-ordinates (because it is a cylinder!).