%%%%%%
% Homoclinc orbits 
% H = 0.5*y^2 - 0.5*x^2 + 0.25*x^4
% (example from the class) 
% dx/dt = dH/dy = y
% dy/dt = -dH/dx = x - x^3
% We are interested in the homoclinc orbit passing by the point H(0,0)=0. 
% It also passes by the point (sqrt(2),0) (H(sqrt(2),0) = 0, and the
% trajectories follow the iso-contours. 

% If your matlab version does not support ode89, replace it by ode45.
% Changing from ode89 to ode45, or changing the total time simulated
% (tspan), or changing the precision (options) will likely change the
% trajectory that is simulated, due to the numerical error : the homoclonic
% orbit is a line in the 2D space, so any numerical error may throw us away
% from it ! 

% Starting position : x = sqrt(2); y = 0; 
y0 = [sqrt(2); 0];

% Max simulated time. You can increase (or decrease it)
max_time = 30;

F = @(t,y) [y(2); y(1)-y(1).^3];
tspan = [0 max_time];
options = odeset('RelTol',1e-13,'AbsTol', 1e-13, 'MaxStep', 1e-3);
[t,y] = ode45(F,tspan,y0, options);

figure(1);
[X,Y] = meshgrid( linspace(-2,2,200),  linspace(-2,2,200));
H = 0.5*Y.^2 - 0.5*X.^2 + 0.25*X.^4;
contourf(X,Y,H, [linspace(-0.25, 0.0, 10) linspace(0.0, 2.0, 10)]); 
title('Contour plot of the Hamiltonian')
grid on; xlabel('x'); ylabel('y'); 
axis equal; 


figure(2);
[X,Y] = meshgrid( linspace(-2,2,200),  linspace(-2,2,200));
H = 0.5*Y.^2 - 0.5*X.^2 + 0.25*X.^4;
contourf(X,Y,H, [linspace(-0.25, 0.0, 10) linspace(0.0, 2.0, 10)]); hold on; 
plot(y(:,1), y(:,2), 'LineWidth', 2, 'Color', 'r')
title('Contour plot of the Hamiltonian, including our simulated trajectory')
grid on; xlabel('x'); ylabel('y'); 
axis equal; 

figure(3); ax = [];
ax(1) = subplot(2,1,1);
plot(t, y(:,1)); hold on; 
grid on
set(gca, 'XTickLabel', []); 
ylabel('Position x')
title('Position (x) and Velocity (y) as a function of time')

ax(2) = subplot(2,1,2);
plot(t, y(:,2)); hold on; 
grid on
xlabel('t (s)'); 
ylabel('Velocity y')
linkaxes(ax, 'x');

% Note : when starting exactly on the homoclinc orbit, the dynamics of your
% system after some time will depend on the numerical error done as your
% approach the fixed point !