%%

% Find one initial point on the local unstable manifold
NumberofPoints=15;%length(LocalUnstablemanifold(:,1))
localunstablemanifold=zeros(NumberofPoints,3);
%options=odeset('e','on','BDF', 'on','RelTol',1e-10,'AbsTol',[1e-10;1e-10;1e-10]);
NumOfBiSection=18;
hx=.0005;

clr=colormap(cool(NumberofPoints+5));

equ=findequilibrium(parameters);
param=parameters;
   xx=equ(1)+hx;
   yy=xx*0+param(3)/(param(6)-1);
    alpha=1-xx-yy;
    beta=1+param(5)*yy./(param(2)+xx);
    zz=(beta.*(a1+alpha)-param(4)*alpha)./(beta-param(4));
    yinit=[yy/5,3*yy/5,yy];
    init=[xx,yy,zz];
        u=yinit(1);
        v=yinit(2);
        w=yinit(3);
        for k=1:NumOfBiSection-1
            init(2)=v;
            [t,solX,a, b, c]=ode15s('chemostatmodel1',[0,-200],init,options,param);
            if c==2
                %yinit(i,1,j)=yinit(i,2,j);
                u=v;
            else
                %yinit(i,3,j)=yinit(i,2,j);
                w=v;
            end
            %yinit(i,2,j)=(yinit(i,1,j)+yinit(i,3,j))/2;
            v=(u+w)/2;
            %yinit(i,4,j)=a;
            %tt(i)=a;
        end
    localunstablemanifold(1,:)=init;
%%
% Find one spiral
   init=localunstablemanifold(1,:);
   init(2)=init(2)+1e-6;
   [~,~,~, b, ~]=ode15s('chemostatmodel3',[0,100], init,options,parameters);

%%
% Find initials of local unstable stable manifold between one spiral return
InitialConditions=localunstablemanifold;

hx1=(b(1)-xx)/NumberofPoints;
    xx=(xx:hx1:(b(1)-hx1));
    yy=xx*0+param(3)/(param(6)-1);
    alpha=1-xx-yy;
    beta=1+param(5)*yy./(param(2)+xx);
    zz=(beta.*(a1+alpha)-param(4)*alpha)./(beta-param(4));
    
    InitialConditions(:,1)=xx'; % along the discrete x-values 
    InitialConditions(:,3)=zz'; % ones(NumberOfRun,1)*equ(3);
    InitialConditions(:,2)=1-InitialConditions(:,1)-InitialConditions(:,3); 

    yinit=zeros(NumberofPoints,4);
    for i=1:NumberofPoints
        yinit(i,1:3)=[InitialConditions(i,2)/5,3*InitialConditions(i,2)/5,InitialConditions(i,2)];
    end

u=zeros(NumberofPoints,1);
v=u;
w=v;
tt=w;

    init=InitialConditions;
    %init(:,2)=yinit(:,2,j);
    yy=yinit;
    init(:,2)=yy(:,2);
        u=yy(:,1);
        v=yy(:,2);
        w=yy(:,3);
        tt=yy(:,4);

    for i=1:NumberofPoints
        %yyinit(1,:)=yinit(i,:,j)
        %Unstablemanifold(i,:,j)=init;
        %yinit(1,:,i)=[InitialConditions(i,2)/5,3*InitialConditions(i,2)/5,InitialConditions(i,2)];
        for k=1:NumOfBiSection-1
            init(i,2)=v(i);
            [t,solX,a, b, c]=ode15s('chemostatmodel1',[0,-200],init(i,:),options,param);
            if c==2
                %yinit(i,1,j)=yinit(i,2,j);
                u(i)=v(i);
            else
                %yinit(i,3,j)=yinit(i,2,j);
                w(i)=v(i);
            end
            %yinit(i,2,j)=(yinit(i,1,j)+yinit(i,3,j))/2;
            v(i)=(u(i)+w(i))/2;
            %yinit(i,4,j)=a;
            tt(i)=a;
        end
        %Unstablemanifold(i,:,j)=init(i,:);   
    end
        localunstablemanifold=init;
        yinit=[u,v,w,tt];
%%
Unstablemanifold=localunstablemanifold;
plot3(Unstablemanifold(:,1),Unstablemanifold(:,3),Unstablemanifold(:,2),'b.','markersize',5)

psection=zeros(length(Unstablemanifold(:,1)),3);
tend=800;
itrt=ones(length(Unstablemanifold(:,1)),1);
for i=1:length(Unstablemanifold(:,1))
    init=Unstablemanifold(i,:);
    [t,solX,a, b, c]=ode15s('chemostatmodelus',[0,tend],init,options,parameters);
                itrt(i)=length(t);
                psection(i,:)=b;
                %tend=16;
                %[t,solX]=ode15s('chemostatmodelus',[0,tend],init,options1,parameters);
                strt=1;%1050;
                t=t(strt:end,1);
                N=0;
                xx=solX(strt:end-N,1);
                yy=solX(strt:end-N,2);
                zz=solX(strt:end-N,3);
                %if mod(i,3)==0
                %plot3(xx,zz,yy,'color',clr(i,:),'linewidth',1,'markersize',5)
                plot3(xx,zz,yy,'b-','linewidth',1,'markersize',5)
                hold on
                %end
end
plk=1;
psection1=psection(plk:end,:);

psectionfamily=psection1;

