function [alphas,omegais]=OrrSommerfeld_BL_temporal(Re);

%close all; clear all;
%clear all
N=300;
%% Base flow
cd ../FalknerSkan
OS_Falkner_Skan % Generates Blasius or BL with non-zero pressure gradient
cd ../OrrSommerfeld


% %% Second- and fourth-order derivatives & boundary conditions
% 
 D2 = D*D; 
 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;
 
 H=10;%delta*max(y);
 
 D=D/H;
 D2=D2/H^2;
 D4=D4/H^4;
 
 D2 = D2(2:N,2:N);
 D4=D4(2:N,2:N); y = y(2:N);
 D = D(2:N,2:N);
 
 
 U=U(2:N);
%% Parameters

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