clear all

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

tran=80;
tend=400;

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


xinit=X0;
[~,X]=ode15s(ModelName, [0 tran], xinit, options, param);

figure(10)
hold off
plot3(X(:,1),X(:,2),X(:,3));

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

X1=X(end,:);% Used as new initial point.

%return
%%

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

xinit=X1;

    [t,X,aa,bb,cc]=ode15s(ModelName, [0 tend], xinit, options, param);
    
trn=1; % Skip the first trn more loops.

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

   aa1=aa(trn:end,:);
    bb1=bb(trn:end,:);
    cc1=cc(trn:end,:);
 
        idx=(bb1(:,2)<0); % Collect only those points with y>0

    aa2=aa1(idx,:);
    bb2=bb1(idx,:);
    cc2=cc1(idx,:);
   
    ppnts=bb2;
    
%%
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