function CyclHenon(p)
% Determination of the periodic orbits of the Henon map
% Cyclist relaxation method
% Map x'=f(x) --> Flow y=C[f(x)-x]
% INPUT p = period
% OUTPUT prime cycles with period p
% Last modified: 12/06/2002 (C. Chandre)
global a b
a=1.4;
b=0.3;

PC=genC(p);
Period(1:p,1)=NaN;
dimP=1;
for ic=1:2^p
    x(1:p)=0;
    C=diag(PC(:,ic));
    options=odeset('AbsTol',1e-13,'RelTol',1e-13,'Events',@divr);
    [tout,xout]=ode45(@modflow,[0 100],x,options,p,C);
    if norm(modflow(tout(end),xout(end,:)',p,C))<=1e-7
        ind1=checkNP(xout(end,:),Period,dimP);
        ind2=checkRP(xout(end,:),p);
        if ind1==1 & ind2==1
            period=sprintf('%+.10f\r',xout(end,:))
            dimP=dimP+1;
            Period(:,dimP)=xout(end,:)';
        end
    end
end

%-----------------------------------
function [value,istermination,direction] = divr(t,x,p,C)
value=norm(x)-1e+2;
istermination=1;
direction=1;

%-----------------------------------
function xp=modflow(t,x,p,C)
global a b
Px(1:p-1)=x(2:p);
Px(p)=x(1);
PIx(2:p)=x(1:p-1);
PIx(1)=x(p);
xp=C*(Px-1+a*x.^2'-b.*PIx)';

%-----------------------------------
function vecC = Irand(p)
for it=1:p
    x=rand(1);
    if x<=1/2
        vecC(it)=1;
    else vecC(it)=-1;
    end
end

%-----------------------------------
function PC = genC(p)
dimPC=1;
PC=zeros(p,2^p);
PC(:,1)=ones(p,1);
while dimPC<=2^p-1
    vecC=Irand(p);
    nPC=zeros(dimPC,1);
    for iv=1:dimPC
        nPC(iv)=norm(vecC'-PC(:,iv));
        if min(nPC)>=1e-5
            dimPC=dimPC+1;
            PC(:,dimPC)=vecC';
        end
    end
end

%-----------------------------------
function ind = checkNP(vec,Period,dimP)
for iv=1:dimP
    nP(iv)=min(abs(vec(1)-Period(:,iv)));
end
if min(nP)<=1e-7
    ind=0;
else ind=1;
end

%------------------------------------
function ind = checkRP(vec,p)
if min(abs(vec(1)-vec(2:p)))<=1e-7
    ind=0;
else ind=1;
end

%
% Copyright (c) 2002. C. Chandre
%

% Built December 5, 2002
% Comments? email chandre@cpt.univ-mrs.fr