clear all

a=.25;
eps=.003;
gamma=5;%8;
c=.3;
param=[a,eps,gamma,c];

%%
options=odeset('e','on','RelTol',1e-10,'AbsTol',1e-10,'BDF', 'on');
options=odeset('e','on','RelTol',1e-12,'AbsTol',1e-12);
Tend=400;
r=.00001;

N=50;
cl=colormap(cool(N));
cl2=colormap(copper(N));

c1 = .2000;
c2 = .5000;
  
g1 =8;
g2 =16;


%%
% Nullcline 
%

figure(1)
hold off

v=-.2:.001:1.05;
w=v.*(v-a).*(1-v);
plot3(v,0*v,w,'k-','linewidth',1)
hold on

ww=v/gamma;
plot3(v,0*v,ww,'r-','linewidth',1)

grid on

xlabel('V','fontsize',16)
ylabel('U = V \prime','fontsize',16)
zlabel('W','fontsize',16)
axis([-.2 1 -.2 .2 -.02 .1])
%view(25,20)
view(0,90)

%%
figure(10)
hold off

v=-.2:.001:1.05;
w=v.*(v-a).*(1-v);
plot3(v,0*v,w,'k-','linewidth',1)
hold on

ww=v/gamma;
plot3(v,0*v,ww,'r-','linewidth',1)

grid on

xlabel('V','fontsize',16)
ylabel('U = V \prime','fontsize',16)
zlabel('W','fontsize',16)
axis([-.2 1 -.2 .2 -.02 .1])
%view(25,20)
view(125,40)

%%
figure(2)
hold off

%figure(3)
%hold off

%%
c=c1;
%gamma=(g1+g2)/2;
    param=[a eps gamma c];
equ =[0,0,0];

    %%%%%%%%  unstable eigenvector 
        
        df=numjacob(@threevctrfld,equ',param);
        [aa, bb]=eig(df);
        cc=sum((real(bb)>0));
        sseq=1:length(aa(1,:));
        unstbleigvctr=real(aa(:,sseq((cc>0))));
        unstbleigvctr=sign(unstbleigvctr(2))*unstbleigvctr;
    init=equ'+r*unstbleigvctr;
    init=init';
    
[t,x]=ode15s('threeDmodelunstmfdsearch',[0,Tend], init, options, param);

figure(1)
plot3(x(:,1),x(:,2),x(:,3),'r-.','linewidth',2)
hold on

figure(2)
plot(t,x(:,1),'r-','linewidth',2)
hold on

%%
c=c2;
    param=[a eps gamma c];
    
    %%%%%%%%  unstable eigenvector 
        df=numjacob(@threevctrfld,equ',param);
        [aa, bb]=eig(df);
        cc=sum((real(bb)>0));
        sseq=1:length(aa(1,:));
        unstbleigvctr=real(aa(:,sseq((cc>0))));
        unstbleigvctr=sign(unstbleigvctr(2))*unstbleigvctr;
    init=equ'+r*unstbleigvctr;
    init=init';
    
[t,x]=ode15s('threeDmodelunstmfdsearch',[0,Tend], init, options, param);

figure(1)
plot3(x(:,1),x(:,2),x(:,3),'k-.','linewidth',2)
hold on

figure(2)
plot(t,x(:,1),'r-','linewidth',2)
hold on

%%
for n=1:N;
    c=(c1+c2)/2;
    %gamma=(g1+g2)/2;
     param=[a eps gamma c];
equ =[0,0,0];
    %%%%%%%%  unstable eigenvector 
        df=numjacob(@threevctrfld,equ',param);
        [aa, bb]=eig(df);
        cc=sum((real(bb)>0));
        sseq=1:length(aa(1,:));
        unstbleigvctr=real(aa(:,sseq((cc>0))));
        unstbleigvctr=sign(unstbleigvctr(2))*unstbleigvctr;
    init=equ'+r*unstbleigvctr;
    init=init';
    
[t,x,~,~,cc1]=ode15s('threeDmodelunstmfdsearch',[0,Tend], init, options, param);

figure(1)
plot3(x(:,1),x(:,2),x(:,3),'-','color',cl(n,:),'linewidth',1)

if n>N/2
figure(10)
plot3(x(:,1),x(:,2),x(:,3),'-','color',cl(n,:),'linewidth',1)
end

figure(2)
plot(t,x(:,1),'-','color',cl(n,:),'linewidth',2)

%%%% next iteration parameter

if cc1(1)==1
    c1=c;
else
    c2=c;
end

end

figure(1)
plot3(x(:,1),x(:,2),x(:,3),'-','color',cl(n,:),'linewidth',2)

%figure(10)
%plot3(x(:,1),x(:,2),x(:,3),'r-','linewidth',3)

usinit=init;

figure(2)
%axis([0,100, -.4 1.2])
xlabel('Time','fontsize',16)
ylabel('V','fontsize',16)

%%
[t,x,aa,bb,cc]=ode15s('threeDmodelunstmfdsearchtwo',[0,Tend], usinit, options, param);

ustend=aa;
psectionpts=bb;

%%
figure(3)
hold off

plot3(v,0*v,w,'k-','linewidth',1)
hold on

ww=v/gamma;
plot3(v,0*v,ww,'r-','linewidth',1)

plot3(x(:,1),x(:,2),x(:,3),'-','color',cl(n,:),'linewidth',2)
%plot3(x(:,1),x(:,2),x(:,3),'k-','linewidth',2)

grid on

xlabel('V','fontsize',16)
ylabel('U = V \prime','fontsize',16)
zlabel('W','fontsize',16)
axis([-.2 1 -.2 .2 -.02 .1])
view(108,45)

figure(10)
plot3(x(:,1),x(:,2),x(:,3),'m-','linewidth',3)

figure(4)
hold off
plot(t,x(:,1),'-','color',cl(n,:),'linewidth',2)
hold on

axis([0,100, -.4 1.2])
xlabel('Time','fontsize',16)
ylabel('V','fontsize',16)

%%
figure(5)
hold off

plot3(v,0*v,w,'k-','linewidth',1)
hold on

ww=v/gamma;
plot3(v,0*v,ww,'r-','linewidth',1)

plot3(x(:,1),x(:,2),x(:,3),'-','color',cl(n,:),'linewidth',2)
%plot3(x(:,1),x(:,2),x(:,3),'k-','linewidth',2)

grid on

xlabel('V','fontsize',16)
ylabel('U = V \prime','fontsize',16)
zlabel('W','fontsize',16)
axis([-.2 1 -.2 .2 -.02 .1])
view(108,45)


figure(6)
hold off
plot(t,x(:,1),'-','color',cl(n,:),'linewidth',2)
hold on

axis([0,100, -.4 1.2])
xlabel('Time','fontsize',16)
ylabel('V','fontsize',16)

%%

rr=.00001;
equ =[0,0,0];
    %%%%%%%%  unstable eigenvector 
        df=numjacob(@threevctrfld,equ',param);
        [aa, bb]=eig(df);
        cc=sum((real(bb)<0));
        sseq=1:length(aa(1,:));
        stbleigvctr=real(aa(:,sseq((cc>0))));
        v1=rr*stbleigvctr(:,2);
        v2=rr*stbleigvctr(:,1);

        s1= -1.1e-4;
        s2=-1e-4;
        s=s1;
        init=1*v1+s*(v2);
        init=init';
            [t,xx,~,~,cc2]=ode15s('threeDmodelstmfdsearch',[0,-Tend], init, options, param);
            plot3(xx(:,1),xx(:,2),xx(:,3),'r-','linewidth',2)
        s=s2;
        init=1*v1+s*(v2);
        init=init';
            [t,xx,~,~,cc2]=ode15s('threeDmodelstmfdsearch',[0,-Tend], init, options, param);
            plot3(xx(:,1),xx(:,2),xx(:,3),'r-','linewidth',2)
        vvals=[];
        uvals=[];
        svals=[];
        KK=30; 
        cl2=colormap(winter(KK));
        Init=zeros(KK,3);
        for kk=1:KK
            s=(s1+s2)/2;
        init=1*v1+s*(v2);
        init=init';
        Init(kk,:)=init;
        %init=v1';
        %plot3([0 init(1)],[0 init(3)],[0 init(3)],'k-','linewidth',1)
 
            [~,~,~,~,cc2]=ode15s('threeDmodelstmfdsearch',[0,-Tend], init, options, param);

            if cc2==1
                s1=s;
            else
                s2=s;
            end
            %options=odeset('e','off','RelTol',1e-8,'AbsTol',1e-8,'BDF', 'on');
            [t,xx,aa2,bb2,cc2]=ode15s('threeDmodelstmfdsearchtwo',[0,-Tend], init, options, param);
            figure(3)
            plot3(xx(:,1),xx(:,2),xx(:,3),'-','color',cl2(kk,:),'linewidth',1)
            figure(4)
plot(ustend+abs(t(end))-abs(t),xx(:,1),'-','color',cl2(kk,:),'linewidth',1)

            if cc2==3
               vvals=[vvals,bb2(1)];
               uvals=[uvals,bb2(2)];
               svals=[svals,s];
            end
        end

%%
[uvals,idx]=sort(uvals);
vvals=vvals(idx);
svals=svals(idx);
idx=(uvals-psectionpts(2)>0);
seq=1:length(svals);
upseq=seq(idx);
s1=svals(upseq(1)-1);
s2=svals(upseq(1));

s1=svals(1);
s2=svals(end);
s1pol=sign(uvals(upseq(1)-1)-psectionpts(2));
s2pol=sign(uvals(upseq(1))-psectionpts(2));
LL=30;
cl3=colormap(copper(LL));
vvals=[];
uvals=[];
svals=[];
for ll=1:LL
    s=(s1+s2)/2;
    init=1*v1+s*(v2);
    init=init';
    [t,xx,~,bb3,cc3]=ode15s('threeDmodelstmfdsearchthree',[0,-Tend], init, options, param);
            if sign(bb3(2)-psectionpts(2))==s1pol
                s1=s;
            else
                s2=s;
            end
            figure(3)
            plot3(xx(:,1),xx(:,2),xx(:,3),'-','color',cl3(ll,:),'linewidth',1)
            figure(4)
%plot(ustend+abs(t(end))-abs(t),xx(:,1),'-','color',cl3(ll,:),'linewidth',2)
            if cc3==1
               vvals=[vvals,bb3(1)];
               uvals=[uvals,bb3(2)];
               svals=[svals,s];
            end
end

figure(3)
 plot3(xx(:,1),xx(:,2),xx(:,3),'-','color',cl(n,:),'linewidth',2)
 
figure(10)
 plot3(xx(:,1),xx(:,2),xx(:,3),'m-','linewidth',3)
 
 
 figure(5)
 plot3(xx(:,1),xx(:,2),xx(:,3),'-','color',cl(n,:),'linewidth',2)
 plot3(equ(1),equ(2),equ(3),'o','markerfacecolor','W','markersize',8)

 view(145,45)
 
 figure(1)
plot3(equ(1),equ(2),equ(3),'o','markerfacecolor','W','markersize',8)
%%
[vmin,idx]=min(abs(vvals-psectionpts(1)));
s=svals(idx);
init=1*v1+s*(v2);
        init=init';
        [t,xx,aa2,bb2,cc2]=ode15s('threeDmodelstmfdsearchtwo',[0,-Tend], init, options, param);
        figure(3)
        plot3(xx(:,1),xx(:,2),xx(:,3),'r-','linewidth',2)
        figure(4)
        plot(ustend+abs(t(end))-abs(t),xx(:,1),'-','color',cl(n,:),'linewidth',2)
        figure(6)
        plot(ustend+abs(t(end))-abs(t),xx(:,1),'-','color',cl(n,:),'linewidth',2)
%plot(ustend+abs(t(end))-abs(t),xx(:,1),'k-','linewidth',2)
 vmin       
umin=uvals(idx)-psectionpts(2)
%%
[a,b]=meshgrid(-.5:.25:1.5,-1.5:.5:.5);
S1=1.25e+4;
S2=.5e+4;
theta=0;
beta=0;
v10=cos(theta)*v1+sin(theta)*v2;
v20=-sin(beta)*v1+cos(beta)*v2;
equ=[0 0 0];
V=equ(1)+S1*a*v10(1)+S2*b*v20(1);
U=equ(2)+S1*a*v10(2)+S2*b*v20(2);
W=equ(3)+S1*a*v10(3)+S2*b*v20(3);
figure(10)
mesh(V,U,W,'edgecolor','g')

S1=.8e+4;
S2=-1.25e+4;
theta=-pi/6;
beta=.4;
v10=cos(theta)*v1+sin(theta)*v2;
v20=-sin(beta)*v1+cos(beta)*v2;
equ=[0 0 0];
V=equ(1)+S1*a*v10(1)+S2*b*v20(1);
U=equ(2)+S1*a*v10(2)+S2*b*v20(2);
W=equ(3)+S1*a*v10(3)+S2*b*v20(3);

figure(3)
mesh(V,U,W,'edgecolor','g')
view(123, 33)
return