% Rayleigh solver

% equation de Rayleigh
    N = 100;
    L=1.1; 
    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;
    
        
    Up=1+tanh(y(2:N));
    U2p=-2*sinh(y(2:N))./cosh(y(2:N)).^3;
    
    
    kk=[0.02:0.02:1.02];
    
    %kk=0.3;
    for p=1:length(kk)
    k=kk(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;
    [taux(p), ip]=max(imag(l));
    freq(p)=real(l(ip));
    %figure;
    %plot(real(l),imag(l),'x');
    %xlabel('real(\omega)');
    %ylabel('imag(\omega)');
    
    end;
    %return;
    figure(3);
    hold all;
    plot(kk,taux(1:length(kk)),'o');
    xlabel('k');
    ylabel('imag(\omega)');
    
    figure(4);
    hold all;
    plot(kk,freq(1:length(kk)),'o');
    xlabel('k');