
% iterating the Henon map
% Period 2-orbits: fig 2.3

clear
lw=2;

henon_map = inline('[a-x.^2+b*y x]','a','b','x','y');

% This initial condition does tend to period-2 (a,b)=(1.28,-0.3)
a=1.28
b=-0.3
x0=-1.5, y0=-2, Npts=50,
% This initial condition does NOT tend to period-2 (a,b)=(1.28,-0.3)
x0=-2.25,y0=-2.25,Npts=4,

% This initial condition does tend to period-2 (a,b)=(1.4,-0.3)
%a=1.4
%b=-0.3
%x0=-1.368, y0=2.5, Npts=50,
% This initial condition does NOT tend to period-2 (a,b)=(1.4,-0.3)
%x0=-1.368-0.001, y0=2.5, Npts=15,

index=[0];

X=x0;
Y=y0;
xo=x0;
yo=y0;


for i=1:Npts
 index=[index i];
 H=henon_map(a,b,xo,yo)
 xn=H(1);
 yn=H(2);
 X=[X xn];
 Y=[Y yn];
 xo=xn;
 yo=yn;
end

figure(1)
clf
subplot(2,1,1)
plot(X,Y,'r+','LineWidth',1)
set(gca,'FontSize',[18]);
axis([-2.5 2.5 -2.5 2.5]); 
hold on
% plot dotted lines for orbits
plot(X,Y,'r:','LineWidth',1)
% plot last 2 pts with different color
plot(X(Npts:Npts+1),Y(Npts:Npts+1),'b+','LineWidth',3,'MarkerSize',10)
set(gca,'FontSize',[18]);
xlabel('x_n'); ylabel('y_n');
title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);

% zoom in arround periodic orbit
subplot(2,1,2)
% plot last 2 pts with different color
plot(X(Npts:Npts+1),Y(Npts:Npts+1),'b+','LineWidth',3)
set(gca,'FontSize',[18]);
axis tight;
ax=axis;
ax=[ax(1)-0.1*(ax(2)-ax(1)) ax(2)+0.1*(ax(2)-ax(1)) ax(3)-0.1*(ax(4)-ax(3)) ax(4)+0.1*(ax(4)-ax(3))]; 
hold on
plot(X,Y,'r+','LineWidth',1)
axis(ax);
plot(X(Npts:Npts+1),Y(Npts:Npts+1),'b+','LineWidth',3,'MarkerSize',10)
xlabel('x_n'); ylabel('y_n');
title(['Zoom around period 2 orbit']);

fprintf('Press any key to continue...');
pause

figure(2)
clf

subplot(3,1,1)
plot(index,X,'r-o','LineWidth',lw)
set(gca,'FontSize',[18]);
xlabel('n');ylabel('x_n');
title(['Henon map: a=',num2str(a),', b=',num2str(b),', (x_0,y_0)=(',num2str(x0),',',num2str(y0),')']);
subplot(3,1,2)
plot(index,Y,'b-o','LineWidth',lw)
set(gca,'FontSize',[18]);
xlabel('n');ylabel('y_n');
subplot(3,1,3)
plot(index,X,'r-','LineWidth',lw)
hold on
plot(index,Y,'b-','LineWidth',lw)
set(gca,'FontSize',[18]);
xlabel('n');ylabel('x_n & y_n');

