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$):

 A: [ EDIT: included both $+$ and $-$ curves, interchanged $k$ and $\omega$ axes as per your image ]
Here is the plot for $d=1$, $c_L = 1.98$, $c_T = 1$, $\omega$ from 0 to 14.  Note that we need $\omega \ge c_L k$ for $q$ to be real, so I took $k$ up to $14/c_L$.  The Maple commands were:

eqs:= eval([tan(pd)/tan(qd) + 4*k^2*pq/(k^2-q^2)^2, 
         tan(pd)/tan(q*d) + (4*k^2*p*q/(k^2-q^2)^2)^(-1)],
        {p=sqrt(omega^2/cl^2-k^2), q=sqrt(omega^2/ct^2 - k^2)});
    eqs:= eval(eqs,{d=1,cl=1.98,ct=1});
    with(plots):
    cols:= [blue,black]:
    display([seq(implicitplot(eqs[i],omega=0..14, k= 0 .. omega/1.98 - .01, grid=[50,50], 
      gridrefine=3, crossingrefine=3, signchange=false, colour=cols[i]),i=1..2)]);


A: The standard methods for numerically solving non-polynomial equations should work to find  $k$ in a given interval for a given value of $\omega$.  In Maple I would use the fsolve command for that.   To plot the solutions given intervals for $\omega$ and $k$ I would use the implicitplot command.
A: 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.
A: There is a description of algorithm in the book Ultrasonic Waves in Solid Media by Joseph L. Rose: 
http://books.google.de/books?id=DEtHDJJ-RS4C&pg=PA110&dq=&redir_esc=y#v=onepage&q&f=false
And here are my two implementations in MATLAB:


*

*Using implicit plot:
close all
clear all
clc
cL = 1.98;
cT = 1;
d = 1;
p = @(k,w)sqrt(w.^2/cL^2-k.^2);
q = @(k,w)sqrt(w.^2/cT^2-k.^2);
symmetric = @(w,k)tan(q(k,w)*d)./q(k,w)+4*k.^2.*p(k,w).*tan(p(k,w)*d)./(q(k,w).^2-k.^2).^2;
asymmetric = @(w,k)q(k,w).*tan(q(k,w)*d)+(q(k,w).^2-k.^2).^2.*tan(p(k,w)*d)./(4*p(k,w).*k.^2);
h1 = ezplot(symmetric,[0 6 0 14]);
hold on
h2 = ezplot(asymmetric,[0 6 0 14]);
set(h1, 'Color', 'b');
set(h2, 'Color', 'k')
legend('symmetric','antisymmetric')
set(gca,'YGrid','on')




*Using fzero function:
% function tries to find all roots of disspersion equations for Lamb waves
close all
clear all
clc
cL = 1.98;
cT = 1;
d = 1;
N = 1000;
omega = linspace(0,6,N);

for idx = N:-1:1
    w = omega(idx);
    p = @(k)sqrt(w.^2/cL^2-k.^2);
    q = @(k)sqrt(w.^2/cT^2-k.^2);
    symmetric = @(k)tan(q(k)*d)./q(k)+4*k.^2.*p(k).*tan(p(k)*d)./(q(k).^2-k.^2).^2;
    asymmetric = @(k)q(k).*tan(q(k)*d)+(q(k).^2-k.^2).^2.*tan(p(k)*d)./(4*p(k).*k.^2);
    try
        lb = 0;
        ub = 14;
        bstep = 0.1;
        tmps = findAllZeros(symmetric,lb,ub,bstep);
        tmpa = findAllZeros(asymmetric,lb,ub,bstep);
        result{idx,1} = [w tmps];
        result{idx,2} = [w tmpa];
    catch ME
        disp(ME)
    end
end

%%
figure
hold on
h1 = plot(NaN,NaN,'b.');
h2 = plot(NaN,NaN,'k.');
hold off
for idx=1:N
    if numel(result{idx,1}) > 1
        x1 = result{idx,1}(1);
        y1 = result{idx,1}(2:end);
        x1 = [x1+0*y1 get(h1,'xdata')];
        y1 = [y1 get(h1,'ydata')];
        set(h1,'xdata',x1,'ydata',y1)        
    end
    if numel(result{idx,2}) > 1        
        x2 = result{idx,2}(1);
        y2 = result{idx,2}(2:end);
        x2 = [x2+0*y2 get(h2,'xdata')];
        y2 = [y2 get(h2,'ydata')];
        set(h2,'xdata',x2,'ydata',y2)
        drawnow       
    end
end
xlim([0 6])
ylim([0 14])
set(gca,'YGrid','on')
xlabel('\omega')
ylabel('k')

Where function findAllZeros:
    function x = findAllZeros(fun,lb,ub,bstep)
    % fun - handle to the function
    % lb - a lower bound for the function.
    % ub - an upper bound for the function.
    % bstep - step for iteration

    x = []; % Initializes x.

    for i=lb:bstep:ub
        if sign(fun(i-bstep))~=sign(fun(i+bstep))
            tmp = fzero(fun, i);
            if isreal(tmp) && abs(fun(tmp))<1 % eliminate complex values and discontinuities
                x = [x tmp]; %#ok<AGROW>
            end
        end
    end

    % Make sure that there are no duplicates.
    x = unique(x);
    DUPE = (diff([x NaN]) < 1e-16) | isnan(x);
    x(DUPE) = [];


Another idea is to plot 3D surface:
    close all
    clear all
    clc
    figure
    cL = 1.98;
    cT = 1;
    d = 1;
    p = @(k,w)sqrt(w.^2/cL^2-k.^2);
    q = @(k,w)sqrt(w.^2/cT^2-k.^2);
    symmetric = @(w,k)tan(q(k,w)*d)./q(k,w)+4*k.^2.*p(k,w).*tan(p(k,w)*d)./(q(k,w).^2-k.^2).^2;
    N = 1000;
    [ww,kk] = meshgrid(linspace(0,6,N),linspace(0,14,N));
    zz = symmetric(ww,kk);
    surf(ww,kk,zz)
    shading interp
    view(2)
    caxis([-5e-10 5e-10])


