% Neutral stability curve of plane Poiseuille flow
  clear all
  close all
  M = 40 ; b = 0 ; 

   revec = linspace(5500,6500,10);
   avec = linspace(.9,1.1,10);
    for iy = 1:length(revec)
     for ix = 1:length(avec)
      Re = revec(iy);
      k = avec(ix);
      [ee,maxreal] = osefun1(Re,k,M);
      lamb(ix,iy)=maxreal;
     end
    end
    [A,RE]=meshgrid(avec,revec);
%    contour(A,RE,lamb',10)%, [0 0], '-k') ; 
  contour(A,RE,lamb', [0 0],'-k')
  axis square
  grid on
  set(gca,'xtick',[.9:.05:1.1])
  set(gca,'ytick',[5500:250:6500])
  xlabel('k','fontsize',14)
  ylabel('Re','fontsize',14)

