% The Henon map for the course in nonlinear dynamics and chaos
% Vered
clear all; close all;
% solve and plot the Henon map, including stable and unstable manifolds:
% x(n+1)= a+b*y(n)-x(n)*x(n)
% y(n+1)=x(n)

% parameter:
a=0.3; b=-1;
%find fixed points:
if (1-b)*(1-b)+4*a<0
    xstarm=0.;
    xstarp=0.;
    print 'no fixed points'
else
    xstarm=0.5*(b-1-sqrt((1-b)*(1-b)+4*a))
    xstarp=0.5*(b-1+sqrt((1-b)*(1-b)+4*a))
    %xstarm is hyperbolic. find its manifolds
    lamdp=-xstarm+sqrt(xstarm*xstarm+b)
    lamdm=-xstarm-sqrt(xstarm*xstarm+b)
end

N=1250;

% initial conditions:
x(1)=xstarp+0.5; y(1)=xstarp;
for n=1:N
  x(n+1)=a+b*y(n)-x(n)*x(n);
  y(n+1)=x(n);
end
%compute stable and unstable manifolds for xstarm;
% Let s be a parameterization of the manifolds, so that 
% xwp(s)=xstarm+sum wp(i)s^i and ywp(s)=xstarm+sum wp(i)(s/lamdp)^i
% require that F(xwp(s),ywp(s))=(xwp(lamdp*s),ywp(lamdp*s)) 
% This leads to the following equations for the corfficients wp(i)
%                           (and equivalentaly for wm(i) with lamdm).
wp(1)=1;
wm(1)=1;
ncof=50;
for i=2:ncof
    sump=0; summ=0;
    for j=1:i-1
        sump = sump+wp(j)*wp(i-j);
        summ = summ+wm(j)*wm(i-j);
    end
    wp(i)=-sump/(2*xstarm+lamdp^i-b*lamdp^(-i));
    wm(i)=-summ/(2*xstarm+lamdm^i-b*lamdm^(-i));
end
%Now use the series form of the manifolds to find them:
xwp(1)=xstarm; ywp(1)=xstarm;
xwm(1)=xstarm; ywm(1)=xstarm;
s=0.1; ds=0.1;
Nman=250;
for n=1:Nman
  xwp(n+1)=xstarm; xwm(n+1)=xstarm;
  ywp(n+1)=xstarm; ywm(n+1)=xstarm;
  for i=1:ncof
      xwp(n+1)=xwp(n+1)+wp(i)*s^i; ywp(n+1)=ywp(n+1)+wp(i)*(s/lamdp)^i;
      xwm(n+1)=xwm(n+1)+wm(i)*s^i; ywm(n+1)=ywm(n+1)+wm(i)*(s/lamdm)^i;
  end 
  s=s+ds;
end

plot(x,y,'.',xwp,ywp,'-',xwm,ywm,'-.','LineWidth',3)
title('Henon map')
xlabel('x')
ylabel('y')