%This function solve the system of non-linear equations (Eigenvalue Problem) for Longitudinal rolls to get the eigenvalue (growth rate)
% Here RaT is the eigenvalue

function Sol=evalfun_NL(n,asq,S,phi,Q,Le,zeta)

%no. of equations
neq=4;

% Order of the matrices
n2=2*n;

% Matrix corresponding to D=d/dz
j=1; mdu=zeros(n2);
while (2*j-1)<=n2
   mdu(1,2*j)=2*j-1;
   j=j+1;
end
for k=2:n2
    j=1;
    while (k+2*j-1)<=n2
         mdu(k,k+2*j-1)=2*(k+2*j-2);
         j=j+1;
    end
end
% Matrix corresponding to D^{2} 
mdsq=mdu*mdu;

% Matrix corresponding to z
mzu=zeros(n2);
for i=1:n2-1
    mzu(i,i+1)=0.5;
    mzu(i+1,i)=0.5;
end
mzu(2,1)=1;
% Matrix for z^{2}
%mzu2=mzu*mzu;

%%% Here, we will replace mzu by (mzu/2+1/2) due to replacing boundary
%%% conditions at z=+-1, Z=(z+1)/2 
%%% Matrices of order n2 

mu=eye(n2); % Identity matrix

%%% Matrix representation for (Q/6)*(3z^2-6z+2)-1
mF=(Q/6)*(3*(mzu/2+mu/2)^2-3*mzu-mu)-mu;

%%% Matrix representation for (D^2-a^2)
md2=(4*mdsq-asq*mu); 

%%% Matrix representation for a^2*cos(phi)
mM2=asq*cos(phi)*mu;

%%% Matrices of order of n 
u=mu(1:n,1:n); s=zeros(n);
F=mF(1:n,1:n);
d2=md2(1:n,1:n);
M2=mM2(1:n,1:n);

%%% Evaluate the Eigenvalue Problem Ax=sigma*Bx %%%
     %%% Formation of matrix 'A'
     
   A=[d2 s s -(S/Le)*M2;
             s d2 -F zeta*u;
                  -F s 2*d2 Q*u;
                       u -(S/(zeta*Le))*M2 (Q/zeta)*u (2/Le)*d2 ];
                   
    %%% Formation of matrix 'B'
    
    B=[s s -M2 s;
            s s s s;
              s -M2 s s ;
                  s s s s];
              
    
  
  %%% Formulation of boundary condition in the matrix form
  for j=1:neq
    A(j*n-1:j*n,:)=zeros(2,neq*n);
    B(j*n-1:j*n,:)=zeros(2,neq*n); 
 end

for j=1:n
    for k=1:neq
        A(k*n-1,(k-1)*n+j)=1;
        A(k*n,(k-1)*n+j)=(-1)^j; 
    end
end
 
lam=eig(A,B); 
p=1;
for j=1:length(lam)
    if real(lam(j))>0
        tes(p)=lam(j); p=p+1;
    end
end
%sort(tes)'
resol=real(tes);

[Sol,par]= min(abs(resol));


%Here we are calling this function to get the maximum value of Ra over zeta (coupling parameter) by using golden section search method

function [sol,par]=gss_NL(n,asq,S,phi,Q,Le)
GS=(3-sqrt(5))/2;
mp(1)=1E-2; mp(2)=0.2; mp(3)=10;
R=evalfun_NL(n,asq,S,phi,Q,Le,mp(2));
fin=0;
    while fin~=1        
        fir=mp(2)-mp(1); sec=mp(3)-mp(2);
        if fir>sec, gap=1; sgap=3;
        else gap=2; sgap=1;
        end
        apoint=((-1)^gap)*GS*(mp(gap+1)-mp(gap))+mp(2);
        Rx=evalfun_NL(n,asq,S,phi,Q,Le,apoint);        
        if Rx<R, mp(2*gap-1)=apoint;
        else mp(sgap)=mp(2); mp(2)=apoint; R=Rx;
        end
        if mp(3)-mp(1)<1E-10, fin=1; end
    end
    Rmax=zeros(1,3);
    for j=1:3
        Rmax(j)=evalfun_NL(n,asq,S,phi,Q,Le,mp(j));
    end
     [sol,ent]=max(Rmax);
      par=mp(ent);
      
%Here we are calling this function to get the minimum value of Ra over wavenumber by using golden section search method 
%and get the critical value of Ra at all parameter values   


function nlin_NL
n=40;
Q=2;
S=-25;
phi=0.523599;
Le=10;


%%% Golden-section search method to minimize over asq %%%
cont=input('Calculate solution? ','s');
if cont=='y'
    a1=input('Enter asq range: ');
    GS=(3-sqrt(5))/2;
    mp(1)=a1(1); mp(2)=(a1(1)+a1(2))/2; mp(3)=a1(2); 
     R=gss_NL(n,mp(2),S,phi,Q,Le);
    fin=0;
    while fin~=1        
        fir=mp(2)-mp(1); sec=mp(3)-mp(2);
        if fir>sec, gap=1; sgap=3;
        else gap=2; sgap=1;
        end
        apoint=((-1)^gap)*GS*(mp(gap+1)-mp(gap))+mp(2);
        Rx=gss_NL(n,apoint,S,phi,Q,Le);
        if Rx>R, mp(2*gap-1)=apoint;
        else mp(sgap)=mp(2); mp(2)=apoint; R=Rx;
        end
        if mp(3)-mp(1)<1E-8, fin=1; end
    end
    for i=1:3
        Rmin(i)=gss_NL(n,mp(i),S,phi,Q,Le);
    end
    [sigma,parmin]=min(Rmin);
    y=[sigma sqrt(mp(parmin)) S (phi*180)/pi Q Le];
    fprintf('R = %12.8f\n min a = %g\nS = %g\nphi = %g\nQ = %g\nLe = %g\n', y) 
end