function [omegais]=OrrSommerfeld_tanh_temporal(N,Re);

%close all; clear all;

%% grid & diff matrices
%N = 180; % degree of highest Chebyshev polynomial
[D,y]=cheb(N); %from Trefethen, Spectral Methods in Matlab
D2 = D*D; D2 = D2(2:N,2:N);
S = diag([0; 1./(1-y(2:N).^2); 0]);
D4 = (diag(1-y.^2)*D^4 - 8*diag(y)*D^3 - 12*D^2)*S;
D4=D4(2:N,2:N); y = y(2:N);
D = D(2:N,2:N);

%% Base flow

H=20; y=y*H; D=D/H; D2=D2/(H^2); D4=D4/(H^4);
U=0.5*(1+tanh(y/2));

alphas=0.1:0.01:1.5;
%Re=1000;
omegais=zeros(size(alphas));


for ii=1:length(alphas)

%% Constructing the eigenvalue problem

alpha=alphas(ii);
%alpha=1.5

    UU=diag(U);
    II=eye(size(D));
    RR=II*((i*alpha*Re)^(-1));

    L1=UU*D2;
    L2=-alpha^2*UU;
    L3=-diag(D2*U);
    L4=-RR*D4;
    L5=+2*RR*alpha^2*D2;
    L6=-(alpha^4)*RR;

    L=L1+L2+L3+L4+L5+L6;

    F1=D2;
    F2=-alpha^2*II;

    F=F1+F2;

    %% Solve the eigenvalue problem

    [V,lambda]=eig(L,F);
    lambda=diag(lambda);

%% fìnd maximum (less stable or most unstable) lambda
    
    [lambdam,I]=max(imag(lambda));
    
    %% plot eigenspectrum
%     figure(ii)
%     plot(real(lambda),imag(lambda),'b.')
%     hold on
%     grid on
%     ylim([-2 0.5])
%     xlim([0 1])
%     figure(ii)
%     plot(real(lambda(I)),imag(lambda(I)),'ro')
%     hold on
    %% calculate growth rate
    
    lambdam;
    omegai=lambdam*alpha;
    
    omegais(ii)=omegai;
    clear omegai alpha
end

%figure
plot(alphas,omegais)
hold on