I think the handiest trick to dealing with these situations is by noting that, if you have a solution $\varphi$ to the Poisson equation $\nabla^2 \varphi = -\rho/\varepsilon_0$, then you can get new ones by applying the cartesian derivatives $\partial/\partial x$ and so on. (This is because $\partial/\partial x$ commutes with the laplacian; the trick fails in curvilinear coordinates where that's no longer the case.)
Many standard solutions can be seen as arising in this way: for example, the point-dipole field is obtained by taking the point-charge field,
$$
\varphi(\mathbf r) = \frac{q}{4\pi\varepsilon_0}\frac{1}{r},
\qquad\text{obeying}\qquad
\nabla^2\varphi = q\delta(\mathbf r)/\varepsilon_0,
$$
and differentiating it with respect to $x$, giving you
$$
\partial_x\varphi(\mathbf r) = -\frac{q}{4\pi\varepsilon_0}\frac{x}{r^3},
\qquad\text{obeying}\qquad
\nabla^2\partial_x\varphi = q\partial_x\delta(\mathbf r)/\varepsilon_0,
$$
where $\partial_x\delta(\mathbf r) = \delta'(x)\delta(y)\delta(z)$ is a point dipole's representation as a distribution-valued charge density. (The higher multipoles can also be obtained in this way, though if you just use cartesian derivatives then you'll need to be careful with how you handle their linear combinations.)
Similarly, the two-displaced-spheres calculation has a natural home in this formalism, where you take as the starting point the field of a uniform ball of charge:
$$
\varphi(\mathbf r) = \begin{cases}\frac{Q}{4\pi\varepsilon_0}\frac{1}{r} & r>a \\ \frac{Q}{4\pi\varepsilon_0} \frac{r^2}{a^3} & r<a\end{cases},
\qquad\text{obeying}\qquad
\nabla^2\varphi = \rho_0\theta(a-r)/\varepsilon_0,
$$
where you now need to differentiate the Heaviside step function to give you
$$
\frac{\partial}{\partial x} \theta(a-r) = \frac{\partial r}{\partial x} \theta'(a-r) = \frac{x}{r} \delta(a-r),
$$
and the internal derivative $\partial_x r^2 = 2x$ gives you a linear potential and therefore a constant electric field. (I'm possibly fudging some constants. Work this stuff through in detail to check it's right.)
Moreover, it should be clear why this derivative trick works for this situation: if you want to model a uniform ball of charge centered at $(\Delta x,0,0)$, you can say that it's the result of a uniform ball of charge at the origin plus the $\sigma\propto\cos(\theta)$ double-shell thing, or you can do a Taylor expansion in $\Delta x$ to get
$$
\varphi_{\Delta x}(\mathbf r) = \varphi_{0}(\mathbf r) + \Delta x \left[ \frac{\partial}{\partial \Delta x}\varphi_{\Delta x}(\mathbf r) \right]_{\Delta x=0} + \mathcal O((\Delta x)^2),
$$
where the derivative along $\Delta x$ is equivalent to a derivative along $x$ (in the same way that active and passive frame transformations are equivalent) and the neglected terms in $\mathcal O((\Delta x)^2)$ are also thrown away in the hand-wavy example.
OK, so that's all nice, but how do you apply it to your cylinder? Well, you start with a uniform cylinder of charge, and you take the derivative of the charge distribution and its solutions with respect to a coordinate orthogonal to the axis. Happy calculating!