%% Parameters

Pr=1;     % Prandtl number


%% Finite-difference operators

N=40;       %number interior points
dz=1/(N+1); % grid size
Z=zeros(N); % zero matrix
I=eye(N);   % identity matrix


% 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;


%mixed at the top
%D4(1,1)=7;
%D4(N,N)=5;

%mixed at the bottom
D4(1,1)=5;
D4(N,N)=7;



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']


for Ra=1600:40:2000

%kk=3;
kk=[0:0.01:10]; % all k-values

% iterate overa ll k-values
for ik=1:length(kk)
    k=kk(ik);
    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(ik)=max(real(eiv)); % store the largest real part of the computed eigenvalue spectrum and keep track of its index
end; 


figure(4);
hold all;
plot(kk,tx,'o');
grid on;
xlabel('k');
ylabel('Re(s)');
title(['Different values of Ra , Pr = ' num2str(Pr) ]);
end;