clear all

options=odeset('event', 'on', 'BDF','on','AbsTol',1e-12,'RelTol',1e-12); % for solve ode, use =[] for default.

xinit=[0   0   0];
param=[.2 .2 5.7];

tend=500;

    [t,X,aa,bb,cc]=ode15s('rossler', [0 tend], xinit, options, param);

%%    
trn=15; % Skip the first trn returning points


%%%%%%% Poincare map on cross-section f_2(x,y,z)=0.

    aa1=aa(trn+1:end,:); % Times of crossing.
    bb1=bb(trn+1:end,:); % Points of crossing.
    %cc1=cc(trn+1:end,:); % Direction of crossing, usually not needed.
 
    idx=(bb1(:,1)>0); % Collect only those points with x>0
        
    aa2=aa1(idx,:);
    bb2=bb1(idx,:);
    %cc2=cc1(idx,:);
  
    ptime=aa2;
    ppnts=bb2; % Poincare return map's orbit  

%%
% Computet the Lyapunov exponent
%
length(ptime) % length of the return map orbit. 

lyexp=lyapunov(ptime,ppnts) % Lyapunov exponent, log base 2.

return
%%
% Optional below:
%
% Plot the attractor
%
figure(1)
hold off
plot3(X(:,1),X(:,2),X(:,3));
hold on;
plot3(ppnts(:,1),ppnts(:,2),ppnts(:,3),'r.')

xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
zlabel('z','fontsize',16)
grid on

%%
%  Poincare return map is approximately 1-dimensional.

k=1; % 1 - x-coordinate, 2 - y-coord, 3 - z-coord
pmin=min(ppnts(:,k));
pmax=max(ppnts(:,k));
pppnts=[ppnts(1:end-1,k) ppnts(2:end,k)];

figure(2)
hold off
plot([pmin pmax],[pmin pmax],'k-')
hold on

plot(pppnts(:,1),pppnts(:,2),'r.')
xlabel('p_n','fontsize',16)
ylabel('p_{n+1}','fontsize',16)
title('Return Map','fontsize',16)

axis([pmin pmax pmin pmax])


return