%figure(100)
%hold off

% Axis
plot3([0 0],[0 1],[0 0],'-','color',[.5 .5 .5],'linewidth',1)
hold on
plot3([0 1],[0 0],[0 0],'-','color',[.5 .5 .5],'linewidth',1)
plot3([0 0],[0 0],[0 1],'-','color',[.5 .5 .5],'linewidth',1)
% Bound of region
plot3([1 0],[0 1],[0 0],'r-','linewidth',1)
plot3([0 0],[0 1],[1 0],'r-','linewidth',1)
plot3([1 0],[0 0],[0 1],'-','color',[.5 .5 .5],'linewidth',1)
% x-nullcline
[a,b]=meshgrid(0:.025:1, 0:.025:1);
x=0*a+a.*b;z=a-a.*b;
alpha=m1*(a2+x).*(1-x-z);
beta=m1*(a2+x);
gamma=a1+1-x-z;
delta=a2+x;
aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
y=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;
cl=0*x;
%mesh(x,z,y,cl);
mesh(x,z,y,'edgecolor', 'g') 
by=max(max(y));

x=0:.01:1;
z=(m2*x./(a2+x)-1).*(a3)/m3;
idx=(z>0);
xx=x(idx);
ax=xx(1)*.88;
z=z-(1-x);
idx=(z>0);
xx=x(idx);
bx=xx(1);

%y-nullcline
yznl=a3/(m3-1);
[a,b]=meshgrid(0:.05:1, 0:.08:1);
x=ax+(bx-ax)*a;y=by*b;
y=yznl+(by-yznl)*b;
z=(m2*x./(a2+x)-1).*(a3+y)/m3;
%mesh(x,z,y,'edgecolor', 'c') % off

[a,b]=meshgrid(0:.05:1, 0:.25:1);
x=ax+(bx-ax)*a;
y=(yznl-.001)*b;
z=(m2*x./(a2+x)-1).*(a3+y)/m3;
%mesh(x,z,y,'edgecolor', 'c') % off

x=0:.01:1;
x=ax+(bx-ax)*x;
y=0*x+a3/(m3-1);
z=(m2*x./(a2+x)-1).*(a3+y)/m3;
% plot3(x,z-.001,y,'b-','LineWidth',1.0,'markersize',5) % off

%z-nullcline
[a,b]=meshgrid(0:1:1, 0:1:1);
x=.88*a;z=.84*b;
y=0*a+a3/(m3-1);
%mesh(x,z,y,'edgecolor', 'b') % off

% Fold of the x-nullcline
zseq=0:.025:1;
seq=0:.001:1;
xfold=zeros(length(zseq),3);
i=1;
for zz=zseq
    xx=(1-zz)*seq;
    alpha=m1*(a2+xx).*(1-xx-zz);
    beta=m1*(a2+xx);
    gamma=a1+1-xx-zz;
    delta=a2+xx;
    aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
    yy=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;
    [bary,idx]=max(yy);
    xfold(i,:)=[xx(idx),yy(idx),zz];
    i=i+1;
end
plot3(xfold(:,1),xfold(:,3),xfold(:,2),'m-','markersize',5)
plot3(0*xfold(:,1),xfold(:,3),xfold(:,2),'m--','markersize',5)
% Find and plot PDLS point
zz=xfold(:,3)';
xx=0*zz;
alpha=m1*(a2+xx).*(1-xx-zz);
beta=m1*(a2+xx);
gamma=a1+1-xx-zz;
delta=a2+xx;
aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
yytran=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;

n=150;%100;
ytop=.8;
ybot=.0;
h=(ytop-ybot)/n;
yy=ybot:h:ytop;
m=2;%50;
nn=m*length(yy);
srch=zeros(1,length(yy));
ypdls=zeros(length(zz),1)';
zpdls=zeros(length(zz),1)';
bary=xfold(:,2);

f = @(y,z) del*z.*(m3*y./(a3+y)-1)./y./(-1-m3*z./(a3+y));
fun = @(y,z) (m1*(1-y-z)./(a1+1-y-z)-1-m2*y/a2)./(-y-m3*z.*y./(a3+y));

j=1;
for z=zz;
    i=1;
    %%%%%%%
        
        [ysolu,zsolu]=rk4ode(bary(j),yy(1),z,nn,f);                 % Solution on x=0 from bary to ymin
    %%%%%%%
    for k=1:length(yy);
        ysolupart=ysolu(1:k*m);
        zsolupart=zsolu(1:k*m);
        %%%%%%%%%
        
        funval = fun(ysolupart,zsolupart);
        q = simpsons(funval,ysolupart(1),ysolupart(end),[]);   % pdls integral to ysolupart(end)
        %%%%%%%%%
        srch(i)=q; 
        i=i+1;
    end
    %[tsting,idx]=min(abs(srch));
    idx=(srch>0);
    seq=1:length(srch);
    aaa=seq(idx);
    idx=aaa(1);
    %ny=yy(idx)*yytran(j);
    zpdls(j)=zsolu(m*idx);
    ypdls(j)=ysolu(m*idx);
    j=j+1;
end
plot3(zz*0,zpdls,ypdls,'k-','markersize',8)

x=0:0.00001:1;
xpdls=0*zpdls;
for j=1:length(zpdls)
    f=m1*(1-x-ypdls(j)-zpdls(j))./(a1+1-x-ypdls(j)-zpdls(j))-1-m2*ypdls(j)./(a2+x);
    [~,idx]=min(abs(f));
    xx=x(idx);
    xpdls(j)=xx(1);
end
%plot3(xpdls+.002,zpdls,ypdls,'k-','linewidth',3,'markersize',8)

% Outline of x-nullcline
xx=0:.035:1;
zz=0*xx;
alpha=m1*(a2+xx).*(1-xx-zz);
beta=m1*(a2+xx);
gamma=a1+1-xx-zz;
delta=a2+xx;
aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
yy=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;
plot3(xx,zz-.001,yy,'r-','linewidth',1)
ytran=yy(1);
plot3([0 0],[0 0],[ytran 1],'r-','linewidth',1)

xx=0.001:.035:1;
zz=1-xx;
alpha=m1*(a2+xx).*(1-xx-zz);
beta=m1*(a2+xx);
gamma=a1+1-xx-zz;
delta=a2+xx;
aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
yy=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;
%plot3(xx,zz,yy,'r-','linewidth',1)

zz=xx;
xx=0*zz;
alpha=m1*(a2+xx).*(1-xx-zz);
beta=m1*(a2+xx);
gamma=a1+1-xx-zz;
delta=a2+xx;
aa=-m2;bb=m2*gamma-delta+beta;cc=gamma.*delta-alpha;
yy=(-bb+sqrt(bb.^2-4*aa.*cc))./aa/2;
plot3(xx,zz,yy,'r-','linewidth',2)
% y-nullcline on x-nullsurface
barx=a2/(m2-1);
xx=barx+(1-barx)*(0:.001:1); 
alpha=m2*xx./(a2+xx)-1;
beta=1-xx-alpha*a3/m3;
gamma=alpha/m3+1;
A=gamma*m2;
B=gamma.*(a2+xx)-(a1+beta)*m2-m1*gamma.*(a2+xx);
C=m1*beta.*(a2+xx)-(a1+beta).*(a2+xx);
yy=(-B-sqrt(B.^2-4*A.*C))./A/2;
zz=alpha.*(a3+yy)/m3;
idx=(yy>0);
xx=xx(idx);
yy=yy(idx);
zz=zz(idx);
plot3(xx,zz-.01,yy,'m-','LineWidth',1.0,'markersize',5)
% z-nullcline on x-nullsurface
xx=(0:.001:1);
yy=xx*0+a3/(m3-1);
alpha=1-xx-yy;
beta=1+m2*yy./(a2+xx);
zz=(beta.*(a1+alpha)-m1*alpha)./(beta-m1);
idx=(zz>0);
plot3(xx(idx),zz(idx),yy(idx),'-','color',[1 .5 0],'LineWidth',1,'markersize',5)

% equilibrium point
equ=findequilibrium(parameters);
plot3(equ(1),equ(3),equ(2),'ko','markersize',5,'markerfacecolor','w')
