function [s1,s2]=sim_cascadesF1eqF2(X1,X2,kap,Tf)
s1={};
s2={};
for k=1:size(X1,1)
%for simulation:
  k1=kap(k,1);
  k2=kap(k,2);
  k3=kap(k,3);
  k4=kap(k,4);
  k5=kap(k,5);
  k6=kap(k,6);
  k7=kap(k,7);
  k8=kap(k,8);
  k9=kap(k,9);
  k10=kap(k,10);
  k11=kap(k,11);
  k12=kap(k,12);
  %system:
  %  E + S0 ---k1-> ES0                 9+1 ---k1-> 5
  %  ES0 ---k2-> E + S0                 5 ---k2-> 9+1
  %  ES0 ---k3-> E + S1                 5 ---k3-> 9+2
  %  F + S1 ---k4-> FS1                 10+2 ---k4-> 6
  %  FS1 -----k5-----> F + S1           6 -----k5-----> 10+2
  %  FS1 -----k6-----> F + S0           6 -----k6-----> 10+1
  %  S1 + P0 ---k7-> S1P0               2+3 ---k7-> 7
  %  S1P0 ---k8-> S1 + P0               7 ---k8-> 2+3
  %  S1P0 ---k9-> S1 + P1               7 ---k9-> 2+4
  %  F + P1 ---k10--> FP1               10+4 ---k10--> 8
  %  FP1 ---k11-> F + P1                8 ---k11-> 10+4
  %  FP1 ---k12-> F + P0                8 ---k12-> 10+3
  %order of the species: S0 S1 P0 P1 S0E S1F P0S1 P1F E F 
  f=@(t,y)[
  -k1*y(1)*y(9)+k2*y(5)+k6*y(6);
  k3*y(5)-k4*y(2)*y(10)+k5*y(6)-k7*y(2)*y(3)+(k8+k9)*y(7);
  -k7*y(2)*y(3)+k8*y(7)+k12*y(8);
  k9*y(7)-k10*y(4)*y(10)+k11*y(8);
  k1*y(1)*y(9)-(k2+k3)*y(5);
  k4*y(2)*y(10)-(k5+k6)*y(6);
  k7*y(2)*y(3)-(k8+k9)*y(7);
  k10*y(4)*y(10)-(k11+k12)*y(8);
  -k1*y(1)*y(9)+(k2+k3)*y(5);
  -k4*y(2)*y(10)+(k5+k6)*y(6)-k10*y(4)*y(10)+(k11+k12)*y(8)
  ];
  %  
  C1=X1(k,5)+X1(k,9);
  C2=X1(k,1)+X1(k,2)+X1(k,5)+X1(k,6)+X1(k,7);
  C3=X1(k,3)+X1(k,4)+X1(k,7)+X1(k,8);
  C4=X1(k,6)+X1(k,8)+X1(k,10);
  x10=[C2,0,C3,0,0,0,0,0,C1,C4]';
  x20=[C2,0,0,C3,0,0,0,0,C1,C4]';
  [t1{k},s1{k}]=ode23s(f,[0 Tf],x10);
  [t2{k},s2{k}]=ode23s(f,[0 Tf],x20); 
end  


%graphics_toolkit ("fltk")
graphics_toolkit ("gnuplot")
v=[1,2,3,4,9,10]; %for plotting S0 S1 P0 P1 E F 

figure (1);
clf ();
hold on
plot(t1{1},s1{1}(:,v(1))','b','LineWidth',1.5)
plot(t1{1},s1{1}(:,v(2))','g','LineWidth',1.5)
plot(t1{1},s1{1}(:,v(3))','r','LineWidth',1.5)
plot(t1{1},s1{1}(:,v(4))','k','LineWidth',1.5)
plot(t1{1},s1{1}(:,v(5))','c','LineWidth',1.5)
plot(t1{1},s1{1}(:,v(6))','m','LineWidth',1.5)
plot(t1{1},X2(1,v(1))*ones(length(t1{1}),1),'b--','LineWidth',1)
plot(t1{1},X2(1,v(2))*ones(length(t1{1}),1),'g--','LineWidth',1)
plot(t1{1},X2(1,v(3))*ones(length(t1{1}),1),'r--','LineWidth',1)
plot(t1{1},X2(1,v(4))*ones(length(t1{1}),1),'k--','LineWidth',1)
plot(t1{1},X2(1,v(5))*ones(length(t1{1}),1),'c--','LineWidth',1)
plot(t1{1},X2(1,v(6))*ones(length(t1{1}),1),'m--','LineWidth',1)
title('x0:S_0=S_{tot},P_0=P_{tot},E=E_{tot},F=F_{tot}, 1st orthant', "fontsize", 14)
axis([0,Tf,0,max(max(s1{1}(:,v)))])
axis square
xlabel('t', "fontsize", 13);
[legh,objh,outh,outm] =legend('S0 ','S1 ','P0 ','P1 ','E ','F ','S0(X2)','S1(X2)','P0(X2)','P1(X2)','E(X2)','F(X2)',-1);
psn=get(legh,'Position');
psn(3:4)=0.70*psn(3:4);
set(legh,'position',psn);
set(objh,'linewidth',2);
print figure1_ss.pdf
hold off

figure (2);
clf ();
hold on
plot(t2{1},s2{1}(:,v(1))','b','LineWidth',1.5)
plot(t2{1},s2{1}(:,v(2))','g','LineWidth',1.5)
plot(t2{1},s2{1}(:,v(3))','r','LineWidth',1.5)
plot(t2{1},s2{1}(:,v(4))','k','LineWidth',1.5)
plot(t2{1},s2{1}(:,v(5))','c','LineWidth',1.5)
plot(t2{1},s2{1}(:,v(6))','m','LineWidth',1.5)
plot(t2{1},X1(1,v(1))*ones(length(t2{1}),1),'b--','LineWidth',1)
plot(t2{1},X1(1,v(2))*ones(length(t2{1}),1),'g--','LineWidth',1)
plot(t2{1},X1(1,v(3))*ones(length(t2{1}),1),'r--','LineWidth',1)
plot(t2{1},X1(1,v(4))*ones(length(t2{1}),1),'k--','LineWidth',1)
plot(t2{1},X1(1,v(5))*ones(length(t2{1}),1),'c--','LineWidth',1)
plot(t2{1},X1(1,v(6))*ones(length(t2{1}),1),'m--','LineWidth',1)
title('x0:S_0=S_{tot},P_1=P_{tot},E=E_{tot},F=F_{tot}, 1st orthant', "fontsize", 14)
axis([0,Tf,0,max(max(s2{1}(:,v)))])
axis square
xlabel('t', "fontsize", 13);
[legh,objh,outh,outm] =legend('S0 ','S1 ','P0 ','P1 ','E ','F ','S0(X1)','S1(X1)','P0(X1)','P1(X1)','E(X1)','F(X1)',-1);
psn=get(legh,'Position');
psn(3:4)=0.70*psn(3:4);
set(legh,'position',psn);
set(objh,'linewidth',2);
print figure2_ss.pdf
hold off
%
%
%
figure (3);
clf ();
hold on
plot(t1{2},s1{2}(:,v(1))','b','LineWidth',1.5)
plot(t1{2},s1{2}(:,v(2))','g','LineWidth',1.5)
plot(t1{2},s1{2}(:,v(3))','r','LineWidth',1.5)
plot(t1{2},s1{2}(:,v(4))','k','LineWidth',1.5)
plot(t1{2},s1{2}(:,v(5))','c','LineWidth',1.5)
plot(t1{2},s1{2}(:,v(6))','m','LineWidth',1.5)
plot(t1{2},X1(2,v(1))*ones(length(t1{2}),1),'b--','LineWidth',1)
plot(t1{2},X1(2,v(2))*ones(length(t1{2}),1),'g--','LineWidth',1)
plot(t1{2},X1(2,v(3))*ones(length(t1{2}),1),'r--','LineWidth',1)
plot(t1{2},X1(2,v(4))*ones(length(t1{2}),1),'k--','LineWidth',1)
plot(t1{2},X1(2,v(5))*ones(length(t1{2}),1),'c--','LineWidth',1)
plot(t1{2},X1(2,v(6))*ones(length(t1{2}),1),'m--','LineWidth',1)
title('x0:S_0=S_{tot},P_0=P_{tot},E=E_{tot},F=F_{tot}, 2nd orthant', "fontsize", 14)
axis([0,Tf,0,max(max(s1{2}(:,v)))])
axis square
xlabel('t', "fontsize", 13);
[legh,objh,outh,outm] =legend('S0 ','S1 ','P0 ','P1 ','E ','F ','S0(X1)','S1(X1)','P0(X1)','P1(X1)','E(X1)','F(X1)',-1);
psn=get(legh,'Position');
psn(3:4)=0.70*psn(3:4);
set(legh,'position',psn);
set(objh,'linewidth',2);
print figure3_ss.pdf
hold off
%
%
figure (4);
clf ();
hold on
plot(t2{2},s2{2}(:,v(1))','b','LineWidth',1.5)
plot(t2{2},s2{2}(:,v(2))','g','LineWidth',1.5)
plot(t2{2},s2{2}(:,v(3))','r','LineWidth',1.5)
plot(t2{2},s2{2}(:,v(4))','k','LineWidth',1.5)
plot(t2{2},s2{2}(:,v(5))','c','LineWidth',1.5)
plot(t2{2},s2{2}(:,v(6))','m','LineWidth',1.5)
plot(t2{2},X2(2,v(1))*ones(length(t2{2}),1),'b--','LineWidth',1)
plot(t2{2},X2(2,v(2))*ones(length(t2{2}),1),'g--','LineWidth',1)
plot(t2{2},X2(2,v(3))*ones(length(t2{2}),1),'r--','LineWidth',1)
plot(t2{2},X2(2,v(4))*ones(length(t2{2}),1),'k--','LineWidth',1)
plot(t2{2},X2(2,v(5))*ones(length(t2{2}),1),'c--','LineWidth',1)
plot(t2{2},X2(2,v(6))*ones(length(t2{2}),1),'m--','LineWidth',1)
title('x0:S_0=S_{tot},P_1=P_{tot},E=E_{tot},F=F_{tot}, 2nd orthant', "fontsize", 14)
axis([0,Tf,0,max(max(s2{2}(:,v)))])
axis square
xlabel('t', "fontsize", 13);
[legh,objh,outh,outm] =legend('S0 ','S1 ','P0 ','P1 ','E ','F ','S0(X2)','S1(X2)','P0(X2)','P1(X2)','E(X2)','F(X2)',-1);
psn=get(legh,'Position');
psn(3:4)=0.70*psn(3:4);
set(legh,'position',psn);
set(objh,'linewidth',2);
print figure4_ss.pdf
hold off
%
