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 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])
