function VFIkedaFig1
% Velocity field aound the hyperbolic fixed
% point of 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);

% fixed point
xfp=0.5327546229407796;
yfp=0.2468967727110331;
zfp=xfp+i*yfp;

xmin=-0.37;
xmax=1.73;
ymin=-2.25;
ymax=0.9;

y0=[0 ymax ymax ymax ymax -0.5 -1.1 ymin ymin];
x0=[xmin xmin 0.4 1 1.4 xmax xmax 1.1 0];

p=1;
C=eye(2);
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
hold on
xlim([xmin xmax]);
ylim([ymin ymax]);
load Ikeda;
plot(x_ik,y_ik,'g.');
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