function [F,G,bd]=surfacewave(e1,u1,e2,u2,e3,u3)
h = 5;

%  u1=2;
%  e1=1.2^2;
%  
%  u3=1.2;
%  e3=1.5^2;
%  
%  u2=-1;
%  e2=-1.8^2;
 
lambda = 1.8;          % mode wavelength
k0 = 2*pi./lambda;
b=k0*h*sqrt(e2*u2-e1*u1);
c=k0*h*sqrt(e2*u2-e3*u3);
 
for k=1:1000
    bd(1)=0;
    bd(k+1)=bd(k)+5*pi/100;
    aa(k)=(u2/u1)*sqrt(b^2+bd(k)^2);
    bb(k)=(u2/u3)*sqrt(c^2+bd(k)^2);
    
    F(k)=-bd(k)*(aa(k)+bb(k))/(bd(k)*bd(k)+aa(k)*bb(k));

       G(k)=tanh(bd(k));
       
end
% hold on 
% plot(bd(1:1000),F(1:1000));
% plot(bd(1:1000),G(1:1000));
% hold off