%% Parameters

Pr=1;     % Prandtl number


%% Finite-difference operators

N=20;       %number interior points
dz=1/(N+1); % grid size
Z=zeros(N); % zero matrix
I=eye(N);   % identity matrix
zz=[0:dz:1]; %z-axis

% Second derivative operator
D2=-2*eye(N);
D2=full(spdiags(ones(N),1,D2));
D2=full(spdiags(ones(N),-1,D2));
D2=D2/dz^2;


% Fourth derivative operator
D4=6*eye(N);
D4=full(spdiags(-4*ones(N),1,D4));
D4=full(spdiags(-4*ones(N),-1,D4));
D4=full(spdiags(ones(N),2,D4));
D4=full(spdiags(ones(N),-2,D4));

%% Boundary conditions

%%rigid
 D4(1,1)=7;
 D4(N,N)=7;

%%stress free
%D4(1,1)=5;
%D4(N,N)=5;

D4=D4/dz^4;

%% Eigenvalue problem:  s M phi = L phi

%    [ D^2-k^2 I ; Z ] [uz']   [Pr(D^2-k^2 I)^2 ; -Pr k^2 I ] [uz']
%  s |               | |   | = |                            | |   |
%    [    Z      ; I ] [th']   [      Ra I      ; D^2-k^2 I ] [th']

%kk=3;
kkk=[0:0.1:10]; % all k-values
Raaa=[0:100:20000];
[kk,Raa]=meshgrid(kkk,Raaa);
%return;
% iterate overa ll k-values
for ik=1:length(kkk)
    ik
    for ia=1:length(Raaa)
    k=kkk(ik);
    Ra=Raaa(ia);
    M=[[D2-k^2*I Z];[Z I ]]; % define L operator
    L=[[Pr*(D4-2*k^2*D2+k^4*I) -Pr*k^2*I];[Ra*I D2-k^2*I]]; % define M operator
    eiv=eig(L,M); % solve the eigenvalue problem
    [tx(ia,ik)]=max(real(eiv)); % store the largest real part of the computed eigenvalue spectrum and keep track of its index
    end; 
end;

figure(5);
contour(kk,Raa,tx);
hold on;
contour(kk,Raa,tx,[0 0],'LineWidth',2);
colorbar

xlabel('k');
ylabel('Ra');