clear all

om=0.48;%0.8,1.2;
e=0.1;
nper=20;
Nt=2001;
T=linspace(0,nper*2*pi,Nt);dt=T(2)-T(1);

        om1=om+e;c1=cos(pi*om1);s1=sin(pi*om1);
        om2=om-e;c2=cos(pi*om2);s2=sin(pi*om2);
        Phi1=[c1 s1/om1;-om1*s1 c1];
        Phi2=[c2 s2/om2;-om2*s2 c2];
        Phi=Phi2*Phi1;
        [u,s]=eig(Phi);
        sd=diag(s);
        [s1,ind]=sort(-abs(sd));
        sd=sd(ind);u=u(:,ind);
        sd1=sd(1);u1=u(:,1);
        sd2=sd(2);u2=u(:,2);
        
        om1=om+e;c1=cos(om1*dt);s1=sin(om1*dt);
        om2=om-e;c2=cos(om2*dt);s2=sin(om2*dt);
        Phi1=[c1 s1/om1;-om1*s1 c1];
        Phi2=[c2 s2/om2;-om2*s2 c2];
        
        
        
        
        xt=zeros(2,Nt);
        xt(:,1)=u1;
        
%         for ia=1:100
%             
%         T2=input('time = ')
%         T1=rem(T2,2*pi);
%         if T1>pi;a=2;
%         else;
%             a=1;
%         end;
%         a
%         end;
      

        for it=2:Nt;
            
            tt=T(it);
            T1=rem(tt,2*pi);
        if T1>pi;
            dPhi=Phi2;
        else;
            dPhi=Phi1;
         end;

          
          xt(:,it)=dPhi*xt(:,it-1);
            
         figure(50);
         plot(xt(1,1:it),xt(2,1:it),'.r')
         xlabel('$x$','Interpreter','Latex','Fontsize',16)
         ylabel('$\dot x$','Interpreter','Latex','Fontsize',16)
         pause(0.1) 
         drawnow 
            
        end;