It's not too difficult to prove from Maxwell's Equations. First of all the equation refers to the special case of waves with harmonic time dependence. Since you can build any form of time dependence from it, using Fourier series or Fourier integrals that's not much of a restriction, and you can use the exponential form of the series or integral which simplifies the equations a lot. Thus we assume that the electromagnetic field is of the following form
\vec{E}(t,\vec{x})=\vec{E}_0(\vec{x}) \exp(-\mathrm{i} \omega t), \quad \vec{B}=\vec{B}_0(\vec{x}) \exp(-\mathrm{i} \omega t).
Here, I use the usual physicisist's convention with a minus sign in the exponent. Usually engineers use the opposite sign convention, but everything is of course equivalent at the end.
Now we plug this ansatz into Maxwell's equations for free fields, i.e., for vanishing charge and current densities. They read
\vec{\nabla} \cdot \vec{E}=0 \; \Rightarrow \; \vec{\nabla} \cdot \vec{E}_0=0, \qquad(1)
\vec{\nabla} \times \vec{E}+\frac{1}{c} \partial_t \vec{B}=0 \; \Rightarrow \; \vec{\nabla} \times \vec{E}_0 -\mathrm{i} k \vec{B}_0=0, \qquad (2)
\vec{\nabla} \cdot \vec{B}=0 \; \Rightarrow \; \vec{\nabla} \cdot \vec{B}_0=0, \qquad (3)
\vec{\nabla} \times \vec{B} - \frac{1}{c} \partial_t \vec{E}=0 \; \Rightarrow \; \vec{\nabla} \times \vec{B}_0 + \mathrm{i} k \vec{E}_0=0. \qquad (4)
I used the abbreviation k=\omega/c. To get an equation for the electric field alone we take the curl of (2) and use (3):
\vec{\nabla} \times (\vec{\nabla} \times \vec{E}_0)-\mathrm{i} k \vec{\nabla} \times \vec{B}_0=\vec{\nabla} \times (\vec{\nabla} \times \vec{E}_0)-k^2 \vec{E}_0=0. \qquad(5)
Now for cartesian (and only cartesian!) coordinates you can write
\vec{\nabla} \times (\vec{\nabla} \times \vec{E}_0)=\vec{\nabla}(\vec{\nabla} \cdot \vec{E}_0)-\Delta \vec{E}_0=-\Delta \vec{E}_0, \qquad(6)
where in the last step we made use of (1).
Plugging this into (5) we find the Helmholtz equation
\Delta \vec{E}_0+k^2 \vec{E}_0=0.
As emphasized above, to use this form of the equation you must use Cartesian coordinates. To write it in orthonormalized curvilinear coordinates, you either must go back to the equation (5) and work out the "double curl operator" or you rewrite everything in cartesian coordinates first. The latter way is often faster, because it already made use of the divergenceless of the electric field.
Let's do this calculation for cylinder coordinates as this is needed for waveguides with cross sections of circular symmetry (e.g., for the coax cable, which is the most important practical example).
We simply need to express the "curvilinear basis vectors" in terms of cartesian coordinates, i.e., use the definitions
\vec{e}_r=\vec{e}_x \cos \varphi+\vec{e}_y \sin \varphi, \quad \vec{e}_{\varphi}=-\vec{e}_x \sin \varphi+\vec{e}_y \cos \varphi,
In cylinder coordinates the third basis vector is simply the cartesian \vec{e}_z.
Now we write
\vec{E}_0=\vec{e}_r E_r + \vec{e}_{\varphi} E_{\varphi} + \vec{e}_z E_z = \vec{e}_x (E_r \cos \varphi - E_{\varphi} \sin \varphi)+\vec{e}_y (E_r \sin \varphi + E_{\varphi} \cos \varphi) + \vec{e}_z E_z.
Now you can apply the Laplace operator in terms of cylinder coordinates, which is a standard formula you can find in any proper textbook about vector calculus. Without further difficulties. For the z component it's of course trivial, because this is already a Cartesian coordinate. This gives your equation, which is written in somewhat simplified form (an also note that you must write E_z everytwhere):
\frac{1}{r} \frac{\partial}{\partial r} \left (r \frac{\partial E_z}{\partial r} \right ) + \frac{1}{r^2} \frac{\partial^2 E_z}{\partial \varphi^2} + \frac{\partial^2 E_z}{\partial z^2} + k^2 E_z=0. \qquad (7)
Now, since the problem is translation invariant in z direction you can further write
\vec{E}_0(\vec{x})=\vec{E}_1(x,y) \exp(\mathrm{i} k_z z)+\vec{E}_2(x,y) \exp(-\mathrm{i} k_z z).
Plugging this ansatz into (7) you see that both E_{1z} and E_{2z} fulfill your equation with
q^2=k^2-k_z^2.
