a three layer model; calculate the linear growth rate

gistfile1.matlab

syms c epsilon k N gamma

A = [-(-c+0.5*(3+epsilon))*(1+k^2)+k*(1+gamma), -c+0.5*(3+epsilon), 0; ...
    N*(-c+0.5*(1+epsilon)), -(-c+0.5*(1+epsilon))*(k^2+2*N)+k*(gamma-N*(0.5-0.5*epsilon)), N*(-c+0.5*(1+epsilon)); ...
    0, -c*N/epsilon, c*(k^2+N/epsilon)+k*(gamma-0.5*N*(1+epsilon)/epsilon)];

det_A = vpa(det(A));

%%
clear
wavenumberSet  = 0.01:0.02:4;
criticalitySet = 0.01:0.02:2;
c1             = zeros(length(criticalitySet),length(wavenumberSet));
c2             = zeros(length(criticalitySet),length(wavenumberSet));
c3             = zeros(length(criticalitySet),length(wavenumberSet));

N = 0.5;
epsilon = 0.1;

for i_wavenumber = 1:length(wavenumberSet)
    for i_criticality = 1:length(criticalitySet)
        
        k  = wavenumberSet(i_wavenumber);
        gamma = criticalitySet(i_criticality);
        
        poly_A(1) = (0.125*(8.*epsilon*k^4 + 8.*epsilon*k^6 + 8.*k^2*N + 8.*epsilon*k^2*N + 8.*k^4*N + 16.*epsilon*k^4*N + 8.*k^2*N^2))/epsilon;
        poly_A(2) = (0.125*(16.*epsilon*gamma*k^3 - 16.*epsilon*k^4 - 8.*epsilon^2*k^4 + 8.*epsilon*k^5 + 24.*epsilon*gamma*k^5 - 16.*epsilon*k^6 - 8.*epsilon^2*k^6 + 8.*gamma*k*N + 8.*epsilon*gamma*k*N - 16.*k^2*N - 24.*epsilon*k^2*N - 8.*epsilon^2*k^2*N + 4.*k^3*N + 8.*epsilon*k^3*N + 4.*epsilon^2*k^3*N + 16.*gamma*k^3*N + 32.*epsilon*gamma*k^3*N - 16.*k^4*N - 40.*epsilon*k^4*N - 16.*epsilon^2*k^4*N - 4.*k^5*N - 8.*epsilon*k^5*N + 4.*epsilon^2*k^5*N + 8.*gamma*k*N^2 - 16.*k^2*N^2 - 8.*epsilon*k^2*N^2 - 12.*k^3*N^2 - 4.*epsilon*k^3*N^2))/epsilon;
        poly_A(3) = (0.125*(8.*epsilon*gamma^2*k^2 - 28.*epsilon*gamma*k^3 - 12.*epsilon^2*gamma*k^3 + 6.*epsilon*k^4 + 8.*epsilon^2*k^4 + 2.*epsilon^3*k^4 + 16.*epsilon*gamma*k^4 + 24.*epsilon*gamma^2*k^4 - 4.*epsilon*k^5 - 4.*epsilon^2*k^5 - 32.*epsilon*gamma*k^5 - 16.*epsilon^2*gamma*k^5 + 6.*epsilon*k^6 + 8.*epsilon^2*k^6 + 2.*epsilon^3*k^6 - 12.*gamma*k*N - 20.*epsilon*gamma*k*N - 8.*epsilon^2*gamma*k*N + 6.*k^2*N + 14.*epsilon*k^2*N + 10.*epsilon^2*k^2*N + 2.*epsilon^3*k^2*N + 4.*gamma*k^2*N + 8.*epsilon*gamma*k^2*N + 4.*epsilon^2*gamma*k^2*N + 8.*gamma^2*k^2*N + 16.*epsilon*gamma^2*k^2*N + 4.*k^3*N + 6.*epsilon*k^3*N - 8.*epsilon^2*k^3*N - 2.*epsilon^3*k^3*N - 16.*gamma*k^3*N - 48.*epsilon*gamma*k^3*N - 24.*epsilon^2*gamma*k^3*N + 2.*k^4*N + 12.*epsilon*k^4*N + 22.*epsilon^2*k^4*N + 4.*epsilon^3*k^4*N - 8.*gamma*k^4*N - 16.*epsilon*gamma*k^4*N + 8.*epsilon^2*gamma*k^4*N + 8.*k^5*N + 18.*epsilon*k^5*N - 2.*epsilon^3*k^5*N + 10.*k*N^2 + 4.*epsilon*k*N^2 + 2.*epsilon^2*k*N^2 - 4.*gamma*k*N^2 - 4.*epsilon*gamma*k*N^2 - 4.*k^2*N^2 + 4.*epsilon*k^2*N^2 - 12.*gamma*k^2*N^2 - 4.*epsilon*gamma*k^2*N^2 + 22.*k^3*N^2 + 20.*epsilon*k^3*N^2 + 6.*epsilon^2*k^3*N^2 + 2.*k^4*N^2 - 2.*epsilon^2*k^4*N^2))/epsilon;
        poly_A(4) = (0.125*(-12.*epsilon*gamma^2*k^2 - 4.*epsilon^2*gamma^2*k^2 + 6.*epsilon*gamma*k^3 + 8.*epsilon^2*gamma*k^3 + 2.*epsilon^3*gamma*k^3 + 8.*epsilon*gamma^2*k^3 + 8.*epsilon*gamma^3*k^3 - 4.*epsilon*gamma*k^4 - 4.*epsilon^2*gamma*k^4 - 16.*epsilon*gamma^2*k^4 - 8.*epsilon^2*gamma^2*k^4 + 6.*epsilon*gamma*k^5 + 8.*epsilon^2*gamma*k^5 + 2.*epsilon^3*gamma*k^5 + 6.*epsilon*gamma*k*N + 8.*epsilon^2*gamma*k*N + 2.*epsilon^3*gamma*k*N + 6.*gamma*k^2*N + 6.*epsilon*gamma*k^2*N - 10.*epsilon^2*gamma*k^2*N - 2.*epsilon^3*gamma*k^2*N - 8.*epsilon*gamma^2*k^2*N - 8.*epsilon^2*gamma^2*k^2*N - 3.*k^3*N - 7.*epsilon*k^3*N - 5.*epsilon^2*k^3*N - 1.*epsilon^3*k^3*N - 4.*gamma*k^3*N + 4.*epsilon*gamma*k^3*N + 20.*epsilon^2*gamma*k^3*N + 4.*epsilon^3*gamma*k^3*N - 4.*gamma^2*k^3*N - 8.*epsilon*gamma^2*k^3*N + 4.*epsilon^2*gamma^2*k^3*N + 2.*k^4*N + 4.*epsilon*k^4*N + 2.*epsilon^2*k^4*N + 8.*gamma*k^4*N + 18.*epsilon*gamma*k^4*N - 2.*epsilon^3*gamma*k^4*N - 3.*k^5*N - 7.*epsilon*k^5*N - 5.*epsilon^2*k^5*N - 1.*epsilon^3*k^5*N - 3.*k*N^2 - 7.*epsilon*k*N^2 - 5.*epsilon^2*k*N^2 - 1.*epsilon^3*k*N^2 + k^2*N^2 + 7.*epsilon*k^2*N^2 + 7.*epsilon^2*k^2*N^2 + epsilon^3*k^2*N^2 + 4.*gamma*k^2*N^2 + 8.*epsilon*gamma*k^2*N^2 + 4.*epsilon^2*gamma*k^2*N^2 - 4.*k^3*N^2 - 14.*epsilon*k^3*N^2 - 12.*epsilon^2*k^3*N^2 - 2.*epsilon^3*k^3*N^2 + 2.*gamma*k^3*N^2 - 2.*epsilon^2*gamma*k^3*N^2 - 3.*k^4*N^2 - 1.*epsilon*k^4*N^2 + 3.*epsilon^2*k^4*N^2 + epsilon^3*k^4*N^2))/epsilon;
        
        waveSpeed = roots(poly_A);
        c1(i_criticality,i_wavenumber) = waveSpeed(1);
        c2(i_criticality,i_wavenumber) = waveSpeed(2);
        c3(i_criticality,i_wavenumber) = waveSpeed(3);
    end
end

%% plot
wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1);
tmp1   = imag(c1.*wavenumber2d);
tmp2   = imag(c2.*wavenumber2d);
tmp3   = imag(c3.*wavenumber2d);
tmpindex = tmp1(:)>0;
growthRate1 = zeros(size(tmp1));
growthRate2 = zeros(size(tmp1));
growthRate1(tmpindex) = tmp1(tmpindex);
tmpindex = tmp2(:)>0;
growthRate1(tmpindex) = tmp2(tmpindex);
tmpindex = tmp2(:)<0;
growthRate2(tmpindex) = tmp2(tmpindex);
tmpindex = tmp3(:)<0;
growthRate2(tmpindex) = tmp3(tmpindex);

figure
contourf(wavenumberSet,criticalitySet,growthRate1,30)
colorbar
set(gca,'YDir','reverse')
xlabel('wavenumber','fontsize',16)
ylabel('1/Sc','fontsize',16)
set(gca,'fontsize',16)
xlim([0 2])

figure
contourf(wavenumberSet,criticalitySet,growthRate2,30)
colorbar
set(gca,'YDir','reverse')
xlabel('wavenumber','fontsize',16)
ylabel('1/Sc','fontsize',16)
set(gca,'fontsize',16)
xlim([0 2])
%ylim([0.2 0.5])
%% calculate the eigen vectors correspond to phase speed (growthRate>0)
wavenumber2d = repmat(wavenumberSet(:)',length(criticalitySet),1);
tmp1   = c1.*wavenumber2d;
tmp2   = c2.*wavenumber2d;
tmpindex = imag(tmp1(:))>0;
mod1 = zeros(size(tmp1));
mod1(tmpindex) = tmp1(tmpindex);
tmpindex = imag(tmp2(:))>0;
mod1(tmpindex) = tmp2(tmpindex);

wavenumberSet  = 0.01:0.02:4;
criticalitySet = 0.01:0.02:2;
N = 0.5;
epsilon = 0.1;

phi = zeros(size(c1,1),size(c1,2),9);
for i_wavenumber = 1:length(wavenumberSet)
    for i_criticality = 1:length(criticalitySet)
        c = mod1(i_criticality,i_wavenumber);
        k  = wavenumberSet(i_wavenumber);
        gamma = criticalitySet(i_criticality);
A = [-(-c+0.5*(3+epsilon))*(1+k^2)+k*(1+gamma), -c+0.5*(3+epsilon), 0; ...
    N*(-c+0.5*(1+epsilon)), -(-c+0.5*(1+epsilon))*(k^2+2*N)+k*(gamma-N*(0.5-0.5*epsilon)), N*(-c+0.5*(1+epsilon)); ...
    0, -c*N/epsilon, c*(k^2+N/epsilon)+k*(gamma-0.5*N*(1+epsilon)/epsilon)];
        [V D] = eig(A);
        phi(i_criticality,i_wavenumber,1:3) = V(1:3)/V(3);
        phi(i_criticality,i_wavenumber,4:6) = V(4:6)/V(6);
        phi(i_criticality,i_wavenumber,7:9) = V(7:9)/V(9);
    end
end

%%
subplot(3,3,1)
imagesc(real(phi(:,:,1)))
subplot(3,3,2)
imagesc(real(phi(:,:,2)))
subplot(3,3,3)
imagesc(real(phi(:,:,3)))
subplot(3,3,4)
imagesc(real(phi(:,:,4)))
subplot(3,3,5)
imagesc(real(phi(:,:,5)))
subplot(3,3,6)
imagesc(real(phi(:,:,6)))
subplot(3,3,7)
imagesc(real(phi(:,:,7)))
subplot(3,3,8)
imagesc(real(phi(:,:,8)))
subplot(3,3,9)
imagesc(real(phi(:,:,9)))