Numerical computation of the Rayleigh-Lamb curves The Rayleigh-Lamb equations:
$$\frac{\tan (pd)}{\tan (qd)}=-\left[\frac{4k^2pq}{\left(k^2-q^2\right)^2}\right]^{\pm 1}$$
(two equations, one with the +1 exponent and the other with the -1 exponent) where
$$p^2=\frac{\omega ^2}{c_L^2}-k^2$$
and
$$q^2=\frac{\omega ^2}{c_T^2}-k^2$$
show up in physical considerations of the elastic oscillations of solid plates. Here, $c_L$, $c_T$ and $d$ are positive constants. These equations determine for each positive value of $\omega$ a discrete set of real "eigenvalues" for $k$. My problem is the numerical computation of these eigenvalues and, in particular, to obtain curves displaying these eigenvalues. What sort of numerical method can I use with this problem? Thanks.
Edit: Using the numerical values $d=1$, $c_L=1.98$, $c_T=1$, the plots should look something like this (black curves correspond to the -1 exponent, blue curves to the +1 exponent; the horizontal axis is $\omega$ and the vertical axis is $k$):


The Rayleigh-Lamb equations:
$$\frac{\tan (pd)}{\tan (qd)}=-\left[\frac{4k^2pq}{\left(k^2-q^2\right)^2}\right]^{\pm 1}$$
are equivalent to the equations (as Robert Israel pointed out in a comment above)
$$\left(k^2-q^2\right)^2\sin  pd \cos  qd+4k^2pq \cos  pd \sin  qd=0$$
when the exponent is +1, and
$$\left(k^2-q^2\right)^2\cos  pd \sin  qd+4k^2pq \sin  pd \cos  qd=0$$
when the exponent is -1. Mathematica had trouble with the plots because $p$ or $q$ became imaginary. The trick is to divide by $p$ or $q$ in convenience. 
Using the numerical values $d=1$, $c_L=1.98$, $c_T=1$, we divide equation for the +1 exponent by $p$ and the equation for the -1 exponent by $q$. Supplying these to the ContourPlot command in Mathematica I obtained the curves

for the +1 exponent, and

for the -1 exponent.
