N = 100;
    L=5; 
    y=[-L:2*L/N:L]';
    dy=2*L/N;
    
    I=eye(N-1);
    dM=-2*ones(N-1,1);
    D2v=diag(dM);
    
    for ip=2:N-2;
        D2v(ip,ip-1)=1;
        D2v(ip,ip+1)=1;
    end;
    D2v(1,2)=1;
    D2v(N-1,N-2)=1;
    
    D2v=D2v/dy^2;
    
    %L=5; %position of the wall
    R=0.9;    
    
    U=1+R*tanh(y);
    Up=U(2:N);
    U2=-2*R*sech(y).^2.*tanh(y);
    U2p=U2(2:N);
    
    
    %ww=[0.02:0.05:1.02];%-0.01i;
    ww=0.4+0.2i;
    for p=1:length(ww)
    k=ww(p)
    A = (diag(Up)*(D2v-k^2*I)-diag(U2p));
    B = D2v-k^2*I;
    [V,ee] = eig(A,B); %generalized eigenvalue problem
    l=diag(ee)*k;
    tx=max(imag(l)); %gaster approximation
    
    
    w=ww(p);
    A0=-w*D2v;
    A1=diag(Up)*(D2v)-diag(U2p);
    A2=w*I;
    A3=-diag(Up);
    [V,ee] = polyeig(A0,A1,A2,A3); %generalized eigenvalue problem
    l=ee;
    [inu, ip]=min(abs(l-(w-i*tx))); %find closest to gaster approximation
    taux(p)=abs(imag((l(ip))));
    freq(p)=real((l(ip)));
    txgas(p)=tx;
    figure(1); hold all;
    plot(real(l),imag(l),'x');
    xlabel('real(k)');
    ylabel('imag(k)');
    %figure;
    %plot(y(2:N),V(:,ip)); %plot eigenfunction
    %xlabel('y');
    %ylabel('\psi');
    %legend(['\omega=' num2str(ee(ix(1)))],['\omega=' num2str(ee(ix(2)))]);
    end;
    return;
    figure(3);
    hold all;
    plot(ww,taux(1:length(ww)),'o');
    hold on;
    plot(ww,txgas,'x');
    xlabel('\omega');
    ylabel('-imag(k)');
    
    figure(4);
    hold all;
    plot(ww,freq(1:length(ww)),'o');
    xlabel('\omega');
    ylabel('real(k)');