function VFIkedaFig2
% Velocity field aound the hyperbolic fixed
% point of f^3 where f is the Ikeda map
% [Chaos] Chapter "Fixed points, and how to get them"
%         Section "Periodic orbits as extremal orbits"
% Last modified: 12/09/2002 (C. Chandre)
clear;
global a
a=0.9*exp(i*0.4);

p=3;
xfp=-0.1352904060;
yfp=-0.3755883015;

zfp=xfp+i*yfp;
xmin=xfp-1e-1;
xmax=xfp+1e-1;
ymin=yfp-1e-2;
ymax=yfp+2e-2;

options=odeset('RelTol',1e-8,'AbsTol',1e-8);

subplot(1,2,1);
C=eye(2);
y0(1:5)=[-0.3638 -0.36377 -0.363791 -0.3637911 -0.363791091];
y0(6:12)=[-0.3662 -0.3661 -0.366176 -0.3661768 -0.36617678 -0.366176788303 -0.366176788302 ];
x0(1:5)=xmin;
x0(6:12)=xmax;
hold on
xlim([xmin xmax]);
ylim([ymin ymax]);
axis square;
plot(xfp,yfp,'o','MarkerSize',8,'LineWidth',3);
for ix=1:length(y0)
    [tout,yout]=ode45(@modflow,[0 10],[x0(ix) y0(ix)],options,p,C);
    plot(yout(:,1),yout(:,2),'LineWidth',1.5);
end

subplot(1,2,2);
C=[1 0;0 -1];
Nx=4;
y0(1:Nx)=linspace(ymin,ymax,Nx);
y0(Nx+1:2*Nx)=y0(1:Nx);
x0(1:Nx)=xmin;
x0(Nx+1:2*Nx)=xmax;
x0(2*Nx+1:3*Nx)=linspace(xmin,xmax,Nx);
x0(3*Nx+1:4*Nx)=x0(2*Nx+1:3*Nx);
y0(2*Nx+1:3*Nx)=ymin;
y0(3*Nx+1:4*Nx)=ymax;
hold on
xlim([xmin xmax]);
ylim([ymin ymax]);
axis square;
plot(xfp,yfp,'o','MarkerSize',8,'LineWidth',3);
for ix=1:length(y0)
    [tout,yout]=ode45(@modflow,[0 10],[x0(ix) y0(ix)],options,p,C);
    plot(yout(:,1),yout(:,2),'LineWidth',1.5);
end

%-----------------------------------
function zik=ikeda(z)
global a
zik=1+a*z*exp(-6*i/(1+abs(z)^2));

%-----------------------------------
function xp=modflow(t,x,p,C)
global a
xt=x(1)+i*x(2);
for it=1:p
    xf=ikeda(xt);
    xt=xf;
end
xp=C*([real(xf);imag(xf)]-x);

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

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