

  %
  clear all
  tp=2*pi;
  Ne=1001;Nom=6001;
  EP=linspace(0,1,Ne);
  OM=linspace(0.4,2.5,Nom);
  lambda=zeros(Ne,Nom);

  for ia=1:length(EP)

      e=EP(ia);

      for ib=1:length(OM)

          om=OM(ib);
          c1=cos(om*(1+e)*tp/2);s1=sin(om*(1+e)*tp/2);
          c2=cos(om*(1-e)*tp/2);s2=sin(om*(1-e)*tp/2);
          PHI1=[c1,s1/(om*(1+e));-om*(1+e)*s1,c1];
          PHI2=[c2,s2/(om*(1-e));-om*(1-e)*s2,c2];
  
          Phi=PHI2*PHI1;
          mu=max(abs(eig(Phi)));
          lambda(ia,ib)=log(mu)/tp;
          % if lambda(ia,ib)<3.e-16
          %     lambda(ia,ib)=NaN;
          % end
      end
  end

  
  figure(45)
  contourf(OM,EP,lambda);colormap jet;colorbar
  xlabel('$\omega$','Interpreter','latex')
  ylabel('$\epsilon$','Interpreter','latex')
  set(gca,'Fontsize',24)



