%% Parameters

Pr=1;     % Prandtl number
Ra=1800; % Rayleigh 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
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; % we analyze one special value of k, so the loop below is 1:1
%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),nimax]=max(real(eiv)); % store the largest real part of the computed eigenvalue spectrum and keep track of its index
end; 

% plot the spectrum

figure(2); 
hold all;
plot(real(eiv),imag(eiv),'bx');
xlabel('Re(s)');
ylabel('Im(s)');
title(['Eigenvalue:  ', 'k = ' num2str(k) ' , Ra = ' num2str(Ra) ' , Pr = ' num2str(Pr)]);
    
%comparaison with theory for free boundary conditions

nt=20;
for in=1:nt
    kn2=k^2+in^2*pi^2;
    s=max(roots([1 (1+Pr)*kn2 Pr*(kn2^2-Ra*k^2/kn2)]));
    s2=min(roots([1 (1+Pr)*kn2 Pr*(kn2^2-Ra*k^2/kn2)]));
    plot(s,0,'ro');
     plot(s2,0,'rd');
end

%get the eigenvector
%ni=3; % an arbitrary eigenvector
ni=nimax; % the most unstable eigenvector
[eigvec, eigva]=eig(L,M); % solve the eigenvalue problem
 eiguz=[0;eigvec(1:N,ni);0]; % remember that u and T have been stacked in the unkwown vector
 eigT=[0;eigvec(N+1:2*N,ni);0];
 
 figure(3);
 plot(zz,eiguz,'bx');
 hold all;
 plot(zz,eigT,'kx');
 xlabel('z');
 ylabel('u_z, T');
 title(['Eigenvector:  ', 'k = ' num2str(k) ' , Ra = ' num2str(Ra) ' , Pr = ' num2str(Pr)]);


%compare with theory
 plot(zz,sin(ni*zz*pi)*max(eiguz),'bo'); %we normalize to compare
 plot(zz,sin(ni*zz*pi)*max(eigT),'ko');
 xlabel('z');
 ylabel('u_z, T');
 title(['Eigenvector:  ', 'k = ' num2str(k) ' , Ra = ' num2str(Ra) ' , Pr = ' num2str(Pr)]);