% tristable_log.m
% for alpha12 = alpha13, alpha21 = alpha23, alpha31 = alpha32
%function tristable_log = tristable_log(a,b,c,d,e,f)
isstable = zeros(5,5,5);
alpha1 = 10.^(1:1);
alpha2 = 10.^(2:2);
alpha3 = 10.^(3:3);
beta1 = 2;
beta2 = 2;
beta3 = 2;
beta4 = 2;
beta5 = 2;
beta6 = 2;
J = zeros(3,3);

% plotting parameters
mksize = 10;
x = 0;
for ii = 1:length(alpha1)
    for jj = 1:length(alpha2)
        for kk = 1:length(alpha3)
        syms u v w;
        f = alpha1(ii)/(1+v^beta1)+alpha1(ii)/(1+w^beta2)-u;
        g = alpha2(jj)/(1+u^beta3)+alpha2(jj)/(1+w^beta4)-v;
        h = alpha3(kk)/(1+u^beta5)+alpha3(kk)/(1+v^beta6)-w;
        [x,y,z] = solve(f,g,h,u,v,w);
        su = size(x);
        eq = 0;
        for ss = 1:su(1)
            uu = double(x(ss));
            vv = double(y(ss));
            ww = double(z(ss));
            if isreal(uu) && isreal(vv) && isreal(ww)
                %compute jacobian
                J(1,1) = -1;
                J(1,2) = -alpha1(ii)*beta1*vv^(beta1-1)/((1+vv^beta1)^2);
                J(1,3) = -alpha1(ii)*beta2*ww^(beta2-1)/((1+ww^beta2)^2);
                J(2,1) = -alpha2(jj)*beta3*uu^(beta3-1)/((1+uu^beta3)^2);
                J(2,2) = -1;
                J(2,3) = -alpha2(jj)*beta4*ww^(beta4-1)/((1+ww^beta4)^2);
                J(3,1) = -alpha3(kk)*beta5*uu^(beta5-1)/((1+uu^beta5)^2);
                J(3,2) = -alpha3(kk)*beta6*vv^(beta6-1)/((1+vv^beta6)^2);
                J(3,3) = -1;
                ev = eig(J);
                if ev(1)<0 && ev(2)<0 && ev(3)<0
                    eq = eq + 1;
                end 
            end    
        end
            if eq >= 3
            isstable(ii,jj,kk) = 1;
            plot3(log10(alpha1(ii)),log10(alpha2(jj)),log10(alpha3(kk)),'b.','markersize',mksize)
            hold on
            plot3(log10(alpha1(ii)),log10(alpha2(kk)),log10(alpha3(jj)),'b.','markersize',mksize)
            hold on
            plot3(log10(alpha1(jj)),log10(alpha2(ii)),log10(alpha3(kk)),'b.','markersize',mksize)
            hold on
            plot3(log10(alpha1(jj)),log10(alpha2(kk)),log10(alpha3(ii)),'b.','markersize',mksize)
            hold on
            plot3(log10(alpha1(kk)),log10(alpha2(ii)),log10(alpha3(jj)),'b.','markersize',mksize)
            hold on
            plot3(log10(alpha1(kk)),log10(alpha2(jj)),log10(alpha3(ii)),'b.','markersize',mksize)
            hold on
            else
            plot3(log10(alpha1(ii)),log10(alpha2(jj)),log10(alpha3(kk)),'r.','markersize',mksize)
            hold on
            plot3(log10(alpha1(ii)),log10(alpha2(kk)),log10(alpha3(jj)),'r.','markersize',mksize)
            hold on
            plot3(log10(alpha1(jj)),log10(alpha2(ii)),log10(alpha3(kk)),'r.','markersize',mksize)
            hold on
            plot3(log10(alpha1(jj)),log10(alpha2(kk)),log10(alpha3(ii)),'r.','markersize',mksize)
            hold on
            plot3(log10(alpha1(kk)),log10(alpha2(ii)),log10(alpha3(jj)),'r.','markersize',mksize)
            hold on
            plot3(log10(alpha1(kk)),log10(alpha2(jj)),log10(alpha3(ii)),'r.','markersize',mksize)
            hold on
            end
        end
    end
end
%tristable_log = isstable();