clear xi yi xs ys
mu=h^2;

% matrice di linearizzazione in [pi ; 0]
A=[1+mu 1; mu 1]
[V,D]=eig(A);
% assicurarsi che la componente y sia positiva
     if V(2,1) < 0
       V(:,1)=-V(:,1);
     end
     if V(2,2) < 0
       V(:,2)=-V(:,2);
     end
 % autovalori ed autovettori
     lam(1)=D(1,1)
     V(:,1)
     lam(2)=D(2,2)
     V(:,2)
% scegliere un punto lungo la separatrice instabile del linearizzato 
eps=1e-3
if lam(1) > lam(2)
     x1=pi + eps*V(1,1);
     y1=eps*V(2,1);
else
     x1=pi + eps*V(1,2);
     y1=eps*V(2,2);
end

% trovare la sua immagine  
y2=y1-mu*sin(x1); 
x2=x1+y2;
% riportarsi all rettangolo 0< x < 2*pi ; -pi < y < pi
x2=x2-(2*pi)*floor(x2/(2*pi));
y2=y2-(2*pi)*floor(y2/(2*pi)); 
if y2 > pi
y2=y2-2*pi;
end
%segmento che unisce (x1,y1) con la sua immagine
np=1000
yi=y1:(y2-y1)/np:y2;
xi=x1:(x2-x1)/np:x2;
nni=size(xi);ni=nni(2);
% scegliere un punto lungo la separatrice stabile del linearizzato 
eps=1e-3
if lam(1) < lam(2)
     x3=pi + eps*V(1,1);
     y3=eps*V(2,1);
else
     x3=pi + eps*V(1,2);
     y3=eps*V(2,2);
end
% per h=1
%x3=pi+eps(1-sqrt(5)/2)
%y3=eps

% trovare la sua immagine inversa  
x4=x3-y3;
y4=y3+mu*sin(x4);
% riportarsi al solito rettangolo
x4=x4-(2*pi)*floor(x4/(2*pi));
y4=y4-(2*pi)*floor(y4/(2*pi));
if y4 > pi
y4=y4-2*pi;
end
%segmento che unisce (x3,y3) con la sua immagine

ys=y3:(y4-y3)/np:y4;
xs=x3:(x4-x3)/np:x4;
nns=size(xs);ns=nns(2);

%inizializzazione della grafica
hold off
plot(xi,yi,'r')
hold on
plot(xs,ys,'g')
w=0.05;
axis([pi-w pi+w -w w])

% iterazione della mappa e della mappa inversa
nit=4;
for j=1:6
for ii=1:nit
% mappa standard in forma vettoriale
yi=yi-mu*sin(xi);yi=yi-(2*pi)*floor(yi./(2*pi));
xi=xi+yi;xi=xi-(2*pi)*floor(xi./(2*pi));
yi(yi>pi)=yi(yi>pi)-2*pi;
% tracciare la separatrice instabile
%plot(xi,yi,'.r')
for i=2:ni
   if abs(xi(i)-xi(i-1)) < .1 & abs(yi(i)-yi(i-1)) <.1
      plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],'r')
   else
%        plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],':r')
        plot([xi(i-1) xi(i)],[yi(i-1) yi(i)],'.r')
   end
end
% inversa della mappa standard in forma vettoriale
xs=xs-ys;xs=xs-(2*pi)*floor(xs./(2*pi));
ys=ys+mu*sin(xs);ys=ys-(2*pi)*floor(ys./(2*pi));
ys(ys>pi)=ys(ys>pi)-2*pi;
%tracciare la separatrice stabile
%plot(xs,ys,'.g')
for i=2:ns
   if abs(xs(i)-xs(i-1)) < .1 & abs(ys(i)-ys(i-1)) <.1
      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],'g')
   else
%      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],':g')
      plot([xs(i-1) xs(i)],[ys(i-1) ys(i)],'.g')
   end
end
ii
pause
end
% scelta finestra
w=input('finestra?')
if w > 3
  axis([0 2*pi -pi pi])
else
  axis([pi-w pi+w -w w])
end
end





