%% compute eigenspectrum of Blasius boundary layer

%close all; clear all;
clear all
N=100

%% Base flow
cd ../FalknerSkan
OS_Falkner_Skan % Generates Blasius or BL with non-zero pressure gradient
cd ../OrrSommerfeld

%% Boundary conditions


%dU=D*U;
%ddU=D2*U;
%ddddU=D4*U;

y = y(2:N); % originalmente seria D2(1:N+1), mas tira-se os extremos.
D = D(2:N,2:N);
D2 =D2(2:N,2:N); 
D4 =D4(2:N,2:N); 

U = U(2:N);
DU = DU(2:N);
D2U = D2U(2:N);
%D4U=D4*U;
%D4U=D4U(2:N);

figure(1)
hold on
plot(DU,y,'r')
plot(D2U,y,'k')
%plot(ddddU,y,'m')

%ddU=ddU(2:N);
%D=D(2:N,2:N); D2=D2(2:N,2:N); D4=D4(2:N,2:N);
%y=y(2:N);
%U=U(2:N);


%% Parameters

alpha=0.2;
Re=5000;



%% Constructing the eigenvalue problem


    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;
    
    L1=UU*D2;
    L2=-alpha^2*UU;
    L3=-diag(D2U);
    L4=-RR*D4;
    L5=+2*RR*alpha^2*D2;
    L6=-(alpha^4)*RR;

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

    %L=L(2:N,2:N);
    
    F1=D2;
    F2=-alpha^2*II;

    F=F1+F2;

    %F=F(2:N,2:N);
    %% 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(2)
     plot(real(lambda),imag(lambda),'k.')
     hold on
     grid on
     ylim([-2 0.5])
     xlim([0 1])
     figure(2)
     plot(real(lambda(I)),imag(lambda(I)),'ro')
     hold on
    %% calculate growth rate
    
    lambdam;
    omegai=lambdam*alpha;