% Orr-Sommerfeld function
  function maxreal = osefun2(Re,k,M)

% N: number of Gauss-Legendre points
% M: number of (1-x^2)^2*Legendre functions 
% The quadratures are exact if N >= M+4 
% If N < M+4, the results are incorrect

  N = M + 4 ;  [x,w,D] = glxdw(N) ;
  D2 = D*D ; D4 = D2*D2 ; ID = eye(N) ;
  a = k ; b = 0 ;
  k2 = a^2 + b^2  ; 
  BF = diag(1-x.^2) ; D2BF = diag(-2*(1+0*x)) ; W = diag(w) ;

% Legendre polynomials generator
  L0 = 1 + 0*x ; L1 = x ;
  LM = zeros(N,M) ;
  LM(:,1) = L0 ; LM(:,2) = L1 ;
  for k = 1:M-2
   icol = k + 1 ;
   LM(:,icol+1) = (2*k+1)*x.*LM(:,icol)/(k+1) - (k/(k+1))*LM(:,icol-1) ;
  end

% PHI-homogeneous basis
  PHI = zeros(N,M) ;  h0 = (1-x.^2).^2 ;
  for icol = 1:M
   PHI(:,icol) = h0.*LM(:,icol);
  end

% A-operator
  A = PHI'*W*(i*a*BF*(D2-k2*ID) - i*a*D2BF - (1/Re)*(D2-k2*ID)^2)*PHI;

% B-operator
  B = -PHI'*W*(D2-k2*ID)*PHI ;

  [EIVC,EIVL] = eig(A,B) ; lamb = diag(EIVL) ;
  maxreal = max(real(lamb));



