clear all;



%Lyapunov 
e=0.5;

EE=[0.05,0.1,0.2,0.3,0.4,0.5,0.75,1,1.25];
EE=linspace(0.01,1.5,101);
Tp=0.4;

TTp=pi*linspace(0.5,2.3,1001);
lambda=zeros(length(TTp),length(EE));

for ie=1:length(EE);
    ie
    e=EE(ie);
for ip=1:length(TTp);
  Tp=TTp(ip);
Nt=201;

tt=linspace(0,Tp,Nt);
Phi=eye(2);
dt=tt(2)-tt(1);

for it=1:Nt;
    t=tt(it);
    om2t=(1+e*sin(2*pi*t/Tp))^2;
    A=[0,1;-om2t,0];
    Phi=expm(A*dt)*Phi;
    
end;

lambda(ip,ie)=log(max(abs(eig(Phi))))/Tp;

end;

% figure(45)
% plot(TTp(1:ip)/(2*pi),lambda(1:ip),'r')
% xlabel('Period');ylabel('\lambda')
% hold on
% drawnow


v=[0.01,0.1,0.2,0.3,0.4,0.5,0.6];

figure(45)
contourf(TTp/(2*pi),EE,lambda',v);colorbar;
xlabel('T/T_0','fontsize',16);
ylabel('epsilon','fontsize',16);

drawnow


end