clear all

e=0.1;
Ne=40;
E=linspace(0,1,Ne);
No=51;
OM1=linspace(1/2-0.2,1/2+0.2,No);
OM2=linspace(1-0.2,1+0.2,No);
OM3=linspace(1.5-0.2,1.5+0.2,2*No);

OM=[OM1,OM2,OM3];

lambda=zeros(No,Ne);

Nt=200;

tt=linspace(0,2*pi,Nt);
delta=tt(2)-tt(1);

for io=1:length(OM)
    om=OM(io);om2=om^2;

    for ie=1:length(E)
        e=E(ie);

Phi=eye(2);

for it=1:Nt-1

    ti=(tt(it)+tt(it+1))/2;
    c=cos(ti);
    A=[0,1;-om2*(1+e*c)^2,0];
    eA=expm(A*delta);
    Phi=eA*Phi;

end


%det(Phi)
%trace(Phi)

%Phi
s=eig(Phi);
sigma=log(max(abs(s)))/(2*pi);

if abs(sigma)<1.e-10
    lambda(io,ie)=NaN;
else
    lambda(io,ie)=log(max(abs(s)))/(2*pi);
end


    end
end


figure(67)
contourf(OM,E,lambda');colorbar;colormap jet
xlabel('$\omega$','interpreter','latex')
ylabel('$\epsilon$','interpreter','latex')
set(gca,'Fontsize',30)











