function [LocalUnstablemanifold, LocalStablemanifold, GlobalUnstablemanifold, GlobalStablemanifold, a2par, slnkidx, shilnikovorbit, searchonz, srcherror, sorbiterror, shooterror]=searchshilnikovorbit(atp, abt, parameters, PIterate, error, NumberOfRun, hx, SBisection, options, options1)
%
%
%

LocalUnstablemanifold=zeros(NumberOfRun,3,PIterate);
LocalStablemanifold=zeros(1,3,PIterate);  
GlobalUnstablemanifold=zeros(NumberOfRun,3,PIterate); 
GlobalStablemanifold=zeros(1,3,PIterate);
srcherror=zeros(PIterate,3);
a2par=zeros(PIterate,1);
slnkidx=zeros(PIterate,1);
shilnikovorbit=zeros(PIterate,4);
hxseq=zeros(PIterate,1);
sorbiterror=zeros(PIterate,2);
searchonz=zeros(SBisection,3,PIterate); 
shooterror=zeros(PIterate,1);

%options1=odeset('e','off','BDF', 'on','RelTol',1e-18,'AbsTol',[1e-18;1e-16;1e-18]);
%options1=odeset('e','off','BDF', 'on','RelTol',1e-10,'AbsTol',[1e-18;1e-16;1e-16]);
%options1=odeset('e','off','BDF', 'on','RelTol',1e-16,'AbsTol',[1e-16;1e-16;1e-16]);

tend=800;
%%
a2seq =[abt, (atp+abt)/2, atp];

a1=parameters(1);
a2=parameters(2);
a3=parameters(3);
m1=parameters(4);
m2=parameters(5);
m3=parameters(6);
eps=parameters(7);
del=parameters(8);

parameters=[a1 a2 a3 m1 m2 m3 eps del];

%%
%tic
j=1;
Stoperr=1;
while j<=PIterate && Stoperr>error
hxseq(j)=hx;
parameters(2)=a2seq(2);
a2=parameters(2);
a2par(j)=a2seq(2);

equ=findequilibrium(parameters);
%%
% Finding the local stable manifold

xintercept=.95; % Tied with cross-section for unstable manifold 
param=[parameters,equ];

        %%%%%%%%  stable eigenvector 
        df=numjacob(@chemovctfld,equ',param);
        [aa, bb]=eig(df);
        cc=sum((real(bb)<0));
        sseq=1:length(aa(1,:));
        stbleigvctr=real(aa(:,sseq((cc>0))));
        unstbleigvctr1=real(aa(:,2));
        unstbleigvctr2=imag(aa(:,2));
        s=(1-(equ(1)/xintercept+equ(2)+equ(3)))/(stbleigvctr(1)/xintercept+stbleigvctr(2)+stbleigvctr(3));
        %%%%%%%%
    ypt=(equ(2)+s*stbleigvctr(2));
    zpt=(equ(3)+s*stbleigvctr(3));
    xpt=(1-ypt-zpt)*xintercept;
 
LocalStablemanifold(1,:,j)=[xpt ypt zpt];

[~,~,~, b, ~]=ode15s('chemostatmodel0',[0,-tend],LocalStablemanifold(1,:,j),options,parameters);

Stablemanifold=b;
GlobalStablemanifold(1,:,j)=b;
%%

% Find one initial point on the local unstable manifold
equy=parameters(3)/(parameters(6)-1);

A=[unstbleigvctr1(1:2,1) unstbleigvctr2(1:2,1)];
B=A^(-1);
s=B*[hx;0];
xx=equ(1)+hx;
yy=equy;
zz=equ(3)+s(1)*unstbleigvctr1(3,1)+s(2)*unstbleigvctr2(3,1);

Unstablemanifold=zeros(NumberOfRun,3);

Unstablemanifold(1,:)=[xx,yy,zz];  

%%
% Find one spiral
    init=Unstablemanifold(1,:);
    init(2)=init(2)+1e-15;
   [~,~,~, b, ~]=ode15s('chemostatmodel3',[0,tend], init,options,parameters);

%%
% Find initials of local unstable stable manifold between one spiral return

hx1=(b(1)-xx)/NumberOfRun;
    xx=(xx:hx1:(b(1)-hx1));
    %yy=xx*0+equ(2);
    zz=xx*0;
    for i=2:length(xx)
        s=B*[xx(i)-equ(1);0];
        zz(i)=equ(3)+s(1)*unstbleigvctr1(3,1)+s(2)*unstbleigvctr2(3,1);
        Unstablemanifold(i,:)=[xx(i) equ(2) zz(i)];        
    end
    LocalUnstablemanifold(:,:,j)=Unstablemanifold;
 %%
% Find global unstable manifold on one cross-section 

psectionfamily=zeros(length(Unstablemanifold(:,1)),3);

parfor i=1:length(Unstablemanifold(:,1))
    init=Unstablemanifold(i,:);
    [~,~,~, b, ~]=ode15s('chemostatmodelus',[0,tend],init,options,parameters);
    if length(b)==0
        b=init;
        b(1)=(1-b(2)-b(3))*xintercept;
    end
    psectionfamily(i,:)=b;
end
GlobalUnstablemanifold(:,:,j)=psectionfamily;
%%
% Find two orbits on two sides of the plane z=equ(3)
% Bisection the x-interval to line up the z=equ(3) plane. 
%

pzz=psectionfamily(:,3)-Stablemanifold(3);
pyy=psectionfamily(:,2)-Stablemanifold(2);
pxx=psectionfamily(:,1)-Stablemanifold(1);
[~,slnk]=min(sqrt(pxx.^2+pyy.^2+pzz.^2));
%pzz=psectionfamily(:,3)-Stablemanifold(3);
%seq=1:length(pzz);
%idx=(pzz>0);
%ss=seq(idx);
%slnk=ss(1);

if slnk<length(pzz) && pzz(slnk)>0
    xinit(1)=Unstablemanifold(slnk,1);
    xinit(3)=Unstablemanifold(slnk+1,1);
end

if slnk==length(pzz) && pzz(slnk)>0
    xinit(1)=Unstablemanifold(slnk,1);
    xinit(3)=Unstablemanifold(slnk,1)+hx;
end

if slnk>1 && pzz(slnk)<0
    xinit(3)=Unstablemanifold(slnk,1);
    xinit(1)=Unstablemanifold(slnk-1,1);
end

if slnk==1 && pzz(slnk)<0
    xinit(3)=Unstablemanifold(slnk,1);
    xinit(1)=Unstablemanifold(slnk,1)-hx;
end

xinit(2)=(xinit(1)+xinit(3))/2;
%%
for kk=1:SBisection
    xx=xinit(2);
        s=B*[xx-equ(1);0];
        zz=equ(3)+s(1)*unstbleigvctr1(3,1)+s(2)*unstbleigvctr2(3,1);
    init(1)=xx; 
    init(3)=zz; 
    init(2)=equ(2); 
    
    % Forward (global) unstable manifold
    [~,~,~, b, ~]=ode15s('chemostatmodelus',[0,tend],init,options,parameters);

    if length(b)==0
        b=init;
        b(1)=(1-b(2)-b(3))*xintercept;
    end
    pz=b(3)-Stablemanifold(3);

    searchonz(kk,:,j)=b;
 
    if pz>0
        xinit(1)=init(1);
    else
        xinit(3)=init(1);
    end
xinit(2)=(xinit(1)+xinit(3))/2;
end

%%
    xx=xinit(2);
        s=B*[xx-equ(1);0];
        zz=equ(3)+s(1)*unstbleigvctr1(3,1)+s(2)*unstbleigvctr2(3,1);
        
Unstablemanifold(slnk,:)=[xx equy zz];
init=Unstablemanifold(slnk,:);
    [~,~,a, b, ~]=ode15s('chemostatmodelus',[0,tend],init,options,parameters);
    if length(b)==0
        b=init;
        b(1)=(1-b(2)-b(3))*xintercept;
        a=0;
    end
                psectionfamily(slnk,:)=b;
strter=sqrt((xx-equ(1)).^2+(equy-equ(2)).^2+(zz-equ(3)).^2); 

%%
pzz=psectionfamily(:,3)-Stablemanifold(3);
pyy=psectionfamily(:,2)-Stablemanifold(2);
%pxx=psectionfamily(:,1)-Stablemanifold(1);

[aa,idx]=sort(pzz);
bb=pyy(idx);
%cc=pxx(idx);
seq=1:length(aa);
aa=aa+1e-16*seq';

a2test=interp1(aa,bb,0,'linear');
shooterror(j)=a2test;

if a2test>0
    atp=a2;
else
    abt=a2;
end
a2seq =[abt, (atp+abt)/2, atp];

srcherror(j,:)=b-Stablemanifold;
Stoperr=max(abs(b-Stablemanifold));
shilnikovorbit(j,2:end)=Unstablemanifold(slnk,:);
slnkidx(j)=slnk;

%%
%[t,solX]=ode15s('chemostatmodelus',[0,40],Unstablemanifold(slnk,:),options1,parameters);
[t,solX]=ode15s('chemostatmodelus',[0,10],b,options1,parameters);
                 
                xx=solX(:,1);
                yy=solX(:,2);
                zz=solX(:,3);
%d=2000;
%strter=min(sqrt((xx(1:d)-equ(1)).^2+(yy(1:d)-equ(2)).^2+(zz(1:d)-equ(3)).^2)); 
[ender,idx]=min(sqrt((xx-equ(1)).^2+(yy-equ(2)).^2+(zz-equ(3)).^2)); 
sorbiterror(j,:)=[strter, ender];
%shilnikovorbit(j,1)=t(d+idx);
shilnikovorbit(j,1)=a+t(idx);
j=j+1;
end
LocalUnstablemanifold(:,:,j:end)=[];
LocalStablemanifold(:,:,j:end)=[];  
GlobalUnstablemanifold(:,:,j:end)=[]; 
GlobalStablemanifold(:,:,j:end)=[];
srcherror(j:end,:)=[];
a2par(j:end,:)=[];
slnkidx(j:end,:)=[];
shilnikovorbit(j:end,:)=[];
hxseq(j:end,:)=[];
sorbiterror(j:end,:)=[];
searchonz(:,:,j:end)=[];
shooterror(j:end,:)=[];
