clear all
close all
%
% Case: F1=F2=F
%
%Order of the species: S0 S1 P0 P1 S0E S1F P0S1 P1F E F          
%
%The transpose of the stoichiometric matrix:
M=[
 1  0  0  0  -1   0   0    0  1  0;
 0  1  0  0   0  -1   0    0  0  1;
 0  1  1  0   0   0  -1    0  0  0;
 0  0  0  1   0   0   0   -1  0  1;
 0  1  0  0  -1   0   0    0  1  0;
 0  1  0  1   0   0  -1    0  0  0
];
%
%The transpose of the "binomials" matrix:
B=[
 1  0  0  0  -1   0   0    0  1  0;
 0  1  0  0   0  -1   0    0  0  1;
 0  1  1  0   0   0  -1    0  0  0;
 0  0  0  1   0   0   0   -1  0  1;
-1  1  0  0   0   0   0    0 -1  1;
 0 -1 -1  1   0   0   0    0  0  1
];
%
%The matrix for the orthogonal subspace of the stoichiometric subspace:
Mperp=[
1 1 0 0 1 1 1 0 0 0;
0 0 1 1 0 0 1 1 0 0;
0 0 0 0 1 0 0 0 1 0;
0 0 0 0 0 1 0 1 0 1
];
%
%The matrix for the orthogonal subspace of the "binomial" subspace:
Bperp=[
1 1 0 1 1 1 1 1 0 0;
0 0 1 1 0 0 1 1 0 0;
0 1 0 1 1 1 1 1 1 0;
1 0 1 0 1 1 1 1 0 1
];

tic

mixed_cond=signs_sigma(Bperp,Mperp); %We check if Sigma is mixed

%
%If Sigma is mixed, we try to find an orthant where the stoichiometric and the orthogonal of the "binomial" subspaces meet:
%
if mixed_cond==1;
  disp('Sigma is mixed')
  [circuitsBperp,circuitsM]=circuits(Bperp,M);%We start by finding the circuits of Bperp and M
  diff_orthants=orthants(circuitsBperp,circuitsM);%We then find orthants where Bperp and M meet, by inspecting supports and signs of the circuits of both matrices
  [X1,X2,kap]=steady_states(diff_orthants,circuitsBperp,circuitsM,Bperp,M); %We now find different steady states for each orthant
  Tf=70;
  sim_cascadesF1eqF2(X1,X2,kap,Tf);%We simulate each system, to check the different steady states 
else
  disp('Sigma is unmixed')
end

toc
