% Henon map: x_n+1 = 1 + y_n - alpha*(x_n)^2 ; y_n+1 = beta*x_n

clear;

% Declare values for alpha and beta:

alpha = .08;
beta = .3;

% Get fixed points:

x1 = 1/(2*alpha)*((beta-1) + sqrt((1-beta)^2 + 4*alpha));
y1 = beta*x1;

x2 = 1/(2*alpha)*((beta-1) - sqrt((1-beta)^2 + 4*alpha));
y2 = beta*x2;

% Get eigenvalues and eigenvectors of Jacobian matrix DF for
% both fixed points.

DF1 = [-2*alpha*x1 1; beta 0];
DF2 = [-2*alpha*x2 1; beta 0];
[v1,D1] = eig(DF1);
[v2,D2] = eig(DF2);

% Now, plot the eigenvectors as an approximation of the subspaces
% of the nonlinear system.  We know from our analysis that
% D1 has only stable eigenvalues and D2 has one unstable and one
% stable, with the unstable eigenvalue first.  We can thus color
% code the stable and unstable eigenspaces (stable: blue, unstable: red).

xnew1 = [x1 x1 x1] + [-.5 0 .5];
xnew2 = [x2 x2 x2] + [-.5 0 .5];

ynew11 = v1(2,1)/v1(1,1)*(xnew1 - x1) + y1;
ynew12 = v1(2,2)/v1(1,2)*(xnew1 - x1) + y1;
plot(xnew1,ynew11,'b');
hold on;
grid on;
plot(xnew1,ynew12,'b');

ynew21 = v2(2,1)/v2(1,1)*(xnew2 - x2) + y2;
ynew22 = v2(2,2)/v2(1,2)*(xnew2 - x2) + y2;
plot(xnew2,ynew21,'r');
plot(xnew2,ynew22,'b');

% Now, plot stable and unstable manifolds.  Choose initial conditions
% close to each fixed point and iterate the map to find the unstable
% manifold.  To find the stable manifold, you need to iterate the inverse
% of the map (x_n = y_n+1/beta, y_n = x_n+1 - 1 + alpha*(y_n+1/beta)^2),
% which is the same as simulating the system backwards in time (starting
% with y_n and getting y_n-1, y_n-2 etc.).  Try to start stable manifold
% calculations near the stable eigenvector and unstable manifold calculations
% near the unstable eigenvectors.  Note that if the eigenvalue corresponding to
% the eigenvector is negative, then you will be able to compute the manifold
% with only one pass, since it is orientation reversing.  If not, then you
% need to compute one half with a point to one side of the fixed point and
% the other half with a point to the other side of the fixed point.

% Start with fixed point x1,y1 which has (for these values of alpha and beta)
% a negative and positive eigenvalue.  Compute the stable manifold for the 
% negative eigenvalue and eigenvector first.

x1start1 = x1 - .01*v1(1,1);
y1start1 = y1 - .01*v1(2,1);
x1start2 = x1 + .01*v1(1,1);
y1start2 = y1 + .01*v1(2,1);
for n = 2:12
   if (n == 2)
       y_np11 = y1start1;
       x_np11 = x1start1;
       y_np12 = y1start2;
       x_np12 = x1start2;
       stable_manifold11(1,:) = [x1start1 y1start1];
       stable_manifold12(1,:) = [x1start2 y1start2];
   else
       y_np11 = yn1;
       x_np11 = xn1;
       y_np12 = yn2;
       x_np12 = xn2;
   end
    xn1 = y_np11/beta;
    yn1 = x_np11 - 1 + alpha*(y_np11/beta)^2;
    xn2 = y_np12/beta;
    yn2 = x_np12 - 1 + alpha*(y_np12/beta)^2;
    stable_manifold11(n,:) = [xn1 yn1];
    stable_manifold12(n,:) = [xn2 yn2];
end
plot(stable_manifold11(:,1),stable_manifold11(:,2),'b.');
plot(stable_manifold12(:,1),stable_manifold12(:,2),'b.');

% Now, work on the other stable manifold (both sides)

x1start1 = x1 - .01*v1(1,2);
y1start1 = y1 - .01*v1(1,2);
x1start2 = x1 + .01*v1(1,2);
y1start2 = y1 + .01*v1(2,2);

for n = 2:7
   if (n == 2)
       y_np11 = y1start1;
       x_np11 = x1start1;
       y_np12 = y1start2;
       x_np12 = x1start2;
       stable_manifold21(1,:) = [x1start1 y1start1];
       stable_manifold22(1,:) = [x1start2 y1start2];
   else
       y_np11 = yn1;
       x_np11 = xn1;
       y_np12 = yn2;
       x_np12 = xn2;
   end
    xn1 = y_np11/beta;
    yn1 = x_np11 - 1 + alpha*(y_np11/beta)^2;
    xn2 = y_np12/beta;
    yn2 = x_np12 - 1 + alpha*(y_np12/beta)^2;
    stable_manifold21(n,:) = [xn1 yn1];
    stable_manifold22(n,:) = [xn2 yn2];
end
plot(stable_manifold21(:,1),stable_manifold21(:,2),'b.');
plot(stable_manifold22(:,1),stable_manifold22(:,2),'b.');

% And now work on the other fixed point - stable manifold first, then
% work on the unstable manifold.

x2start1 = x2 - .01*v2(1,2);
y2start1 = y2 - .01*v2(1,2);
x2start2 = x2 + .01*v2(1,2);
y2start2 = y2 + .01*v2(2,2);

for n = 2:5
   if (n == 2)
       y_np11 = y2start1;
       x_np11 = x2start1;
       y_np12 = y2start2;
       x_np12 = x2start2;
       stable_manifold31(1,:) = [x2start1 y2start1];
       stable_manifold32(1,:) = [x2start2 y2start2];
else
       y_np11 = yn1;
       x_np11 = xn1;
       y_np12 = yn2;
       x_np12 = xn2;
   end
    xn1 = y_np11/beta;
    yn1 = x_np11 - 1 + alpha*(y_np11/beta)^2;
    xn2 = y_np12/beta;
    yn2 = x_np12 - 1 + alpha*(y_np12/beta)^2;
    stable_manifold31(n,:) = [xn1 yn1];
    stable_manifold32(n,:) = [xn2 yn2];
end
plot(stable_manifold31(:,1),stable_manifold31(:,2),'b.');
plot(stable_manifold32(:,1),stable_manifold32(:,2),'b.');

% And finally the unstable manifold, which is easier...

x2start1 = x2 - .01*v2(1,1);
y2start1 = y2 - .01*v2(2,1);
x2start2 = x2 + .01*v2(1,1);
y2start2 = y2 + .01*v2(2,1);

for n = 2:12
   if (n == 2)
       y_n11 = y2start1;
       x_n11 = x2start1;
       y_n12 = y2start2;
       x_n12 = x2start2;
       unstable_manifold11(1,:) = [x2start1 y2start1];
       unstable_manifold12(1,:) = [x2start2 y2start2];
else
       y_n11 = y_np11;
       x_n11 = x_np11;
       y_n12 = y_np12;
       x_n12 = x_np12;
   end
    x_np11 = 1 + y_n11 - alpha*(x_n11)^2;
    y_np11 = beta*x_n11;
    x_np12 = 1 + y_n12 - alpha*(x_n12)^2;
    y_np12 = beta*x_n12;
    unstable_manifold11(n,:) = [x_np11 y_np11];
    unstable_manifold12(n,:) = [x_np12 y_np12];
end
plot(unstable_manifold11(:,1),unstable_manifold11(:,2),'r.');
plot(unstable_manifold12(:,1),unstable_manifold12(:,2),'r.');
axis([-12 2 -9 1]);
% At this point, you could print out the plot to get an idea
% what the stable and unstable manifolds look like.

% If you want, you can uncomment the following lines to see how
% a few orbits look...

% Now, plot a nearby orbit to each fixed point to give an idea what 
% nearby orbits would do.

%x = [1.8 -9 -10.5];
%y = [-.5 -5 -1.5];
%counter = 0;
%for i = 1:3
%        clear orbit;
%        orbit(1,:) = [x(i) y(i)];
%        xn = x(i);
%        yn = y(i);
%        for n = 1:10
%            if (n ~= 1)
%                xn = xnp1;
%                yn = ynp1;
%            end
%            xnp1 = 1 + yn - alpha*(xn)^2;
%            ynp1 = beta*xn;
%            orbit(n+1,:) = [xnp1 ynp1];
%        end
%        plot(orbit(:,1),orbit(:,2),'k.-');
%end