Plasmon modes of cylinder metalic particle Solving Laplace equation gives plasmon modes of spherical metalic
particle radius $R$, plasma frequency $\omega_p$. Famous result is  $l = 0, 1, 2, ...$
$$\omega_l = \sqrt{\frac{l}{2l+1}} \omega_p$$
What is same result of cylinder metalic particle, radius $R$ and height $H$? Is it known?
 A: For a cylinder (and for most geometries for that matter), there is no analytical solutions for the plasmonic resonances.
This is directly linked to the (in)ability to solve Laplace's equation in the given geometry.
The closest shape that approaches a cylinder where some analytical expressions are available would be an ellipsoid, more particularly a spheroid.
Here I'll follow mainly ref 1.
First, as it is the case for the equation you have for the sphere, we assume that the ellipsoid is small compared to the wavelength so that we are in the quasistatic approximation (when talking about plasmonics or plasmon modes, the quasistatic limit is often implied).
The following results are obtain by solving Laplace's equation for the potential in ellipsoidal coordinate. 
The particle has permittivity $\epsilon_1$ and is surrounded by a medium of permittivity $\epsilon_m$. The three semi axis are $a$, $b$, and $c$. We consider that $a>b>c$.
The polarizability of the ellipsoid along one of the three axis is ($n=a$, $b$ or $c$)
\begin{gather}
\alpha_n=4\pi a b c \frac{\epsilon_1-\epsilon_m}{3\epsilon_m+3L_n(\epsilon_1-\epsilon_m)}
\end{gather}
with
\begin{gather}
L_n=\frac{abc}{2}\int_0^\infty\frac{dq}{(n^2+q)f(q)} \,, \ \ \ \ \ \ \ \ \ \ \ L_a+L_b+L_c=1   
\\
\\
f(q)=\sqrt{(q+a^2)(q+b^2)(q+c^2)} \,\,.
\end{gather}
The dipolar resonance occurs when the polarizability diverges, i.e. when its denominator vanishes. If we consider metallic particles, we have a frequency dependent permittivity $\epsilon_1(\omega)$ that will be the free parameter in the polarizability equation.
Thus we have a resonance when
\begin{equation}
\epsilon_1(\omega)=\bigg(1-\frac{1}{L_n}\bigg)\epsilon_m
\end{equation}
For a sphere, we have $L=1/3$ and we retrieve the known equation $\epsilon_1(\omega)=-2\epsilon_m$, which leads to the equation you have when $l=1$ and a Drude model is used for $\epsilon_1$. 
If the ellipsoid is a spheroid, either prolate (b=c) or oblate (a=b), there is two analytical expressions for the $L_a$ coefficient (the on corresponding to the longer axis):
\begin{gather}
L_a^{pr}=\frac{1-e^2}{e^2}\bigg[-1+\frac{1}{2e}\text{ln}\bigg(\frac{1+e}{1-e}\bigg)\bigg] \ \ \ \ \ \ \ \ \ \ \ \ \ \ e^2=1-\frac{b^2}{a^2}
\\
\\
L_a^{ob}=\frac{g(e)}{2e^2}\bigg[\frac{\pi}{2}-\frac{1}{\tan\big(g(e)\big)}\bigg]-\frac{g^2(e)}{2}
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
g(e)=\sqrt{\frac{1-e^2}{e^2}} \ \ \ \ \ \ \ \ e^2=1-\frac{c^2}{a^2}
\end{gather}
The other $L$ can be directly found because of the relation $L_a+L_b+L_c=1$.
Below is a plot of the two resonant energies of gold nano-spheroids in water for a range of long semi axis $a$ values in nm, computed with the above equations; $\omega_1$, respectively $\omega_2$ are the resonances along the long and short axis.

I have found this paper (https://arxiv.org/abs/0811.4070) where they obtain results for higher order modes.
1 Bohren and Huffman "Absorption and Scattering of Light by Small Particles", pp 141-146
