% Homoclinic orbit example:
% dx/dt = y;   dy/dy = x-3ax^2

% The Hamiltonian is
% H(x,y) = (y^2 - x^2)/2+a*x^3

clear

a=1.0;
fhomo = @(t,y) [y(2); y(1)-3*a*(y(1)).^2];

%[t,x]=ode45(fduff,[0:2*pi/om:5000/om],[1,0]);

% Add quiverplot!

[x,y] = meshgrid(-0.5:0.1:0.5,-0.5:0.1:0.5);
dx = y;
dy = x-3*a*x.^2;
figure(1)
% Normalize
r = ( dx.^2 + dy.^2 ).^0.5;
h=quiver(x,y,dx./r,dy./r);
set(h,'AutoScale','on', 'AutoScaleFactor',0.5)
axis([-.5 .5 -.5 .5])

figure(2)
h=quiver(x,y,dx./r,dy./r);
set(h,'AutoScale','on', 'AutoScaleFactor',0.5)
axis([-.5 .5 -.5 .5])

hold on
for k=1:5
    [t,x]=ode45(fhomo,[0:0.1:20],[.1*k,0]);
    plot(x(:,1),x(:,2),'k')
    [t,x]=ode45(fhomo,[0:0.1:2.85],[-.1*k,.5]);
plot(x(:,1),x(:,2),'k')
end

% Separatrices y = +/- x* sqrt(1-2ax) from H(x,y)=0.
x=-.5:.01:.5;
plot(x,x.*sqrt(1-2*a*x),'g')
plot(x,-x.*sqrt(1-2*a*x),'g')

hold off
axis square

