Either do it directly, without Stokes theorem or use Stokes theorem to find the integral over the entire cylinder, then subtract of the integral over the "flat" ends.
To integrate directly over the cylinder [itex]\{(x,y,z)| x^2+ y^2= 9, 0\le z\le 5}[/itex] use cylindrical coordinates with r set equal to 3, the radius of the cylinder. Then [itex]x= 3cos(\theta)[/itex], [itex]y= 3sin(\theta)[/itex], z= z. You can write the "position vector" of a point on that surface as [itex]\vec{r}= 3cos(\theta)\vec{i}+ 3sin(\theta)\vec{j}+ z\vec{k}[/itex]. The derivatives of that are [itex]\vec{r}_\theta= -3sin(\theta)\vec{i}+ 3cos(\theta)\vec{j}[/itex] and [itex]\vec{r}_z= \vec{k}[/itex] and are vectors in the tangent plane to the cylinder at each point. The cross product [itex]\vec{r}_\theta\times\vec{r}_z= 3cos(\theta)\vec{i}+ 3sin(\theta)\vec{j}[/itex] gives the "vector differential of surface area", [itex]3(cos(\theta)\vec{i}+ sin(\theta)\vec{j})d\theta dz[/itex]. Take the dot product of that with your vector function and integrate.
Since here your vector function is [itex]2xy\vec{i}- y^2\vec{j}+ (z+ xy)\vec{k}[/itex][itex]= 2sin(\theta)cos(\theta)\vec{i}- sin^2(\theta)\vec{j}+ (z+ sin(\theta)cos(\theta))\vec{k}[/itex], that dot product is [itex]6cos^2(\theta)sin(\theta)- 3sin^3(\theta) and the integral would be <br />
[tex]3\int_{z= 0}^5\int_{\theta= 0}^{2\pi}(2cos^2(\theta)sin(\theta)- sin^2(\theta))d\theta dz[/tex]<br />
[tex]= 3\int_{z=0}^5\int_{\theta= 0}^{2\pi}(2cos^2(\theta)0- sin^2(\theta))sin(\theta)d\theta dz[/tex]<br />
[tex]= 3\int_{z=0}^5\int_{\theta= 0}^{2\pi}(3cos^2(\theta)- 1)sin(\theta)d\theta dz[/tex][/itex]