#!/usr/bin/octave -q
#RESULTS
#The naive dynamical systems have 5 steady states each, which are (left discrete, right continuous):
#	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]	[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
#	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]	[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99896, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.99896, 0.99997, 0.0, 0.0, 1.0, 0.0, 0.9707, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
#	[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1]	[0.0, 0.0, 0.0, 0.0, 0.71487, 0.97215, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00485, 0.0, 0.0, 0.0, 0.89503, 5.0E-4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.89503]
#	[0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0]	[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.00839, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.93037, 0.99719, 0.00209, 2.6E-4, 0.0, 0.0, 3.0E-5, 0.0, 0.04366, 0.0, 0.0, 0.93037, 0.0, 0.0, 0.99991, 0.0, 0.0, 0.0, 0.0]
#	[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]	[1.6035E-4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

#Nodes have the following order:
node_names = ["Foxp3";"GATA3";"IFNb";"IFNbR";"IFNg";"IFNgR";"IL10";"IL10R";"IL12";"IL12R";"IL17";"IL18";"IL18R";"IL2";"IL23";"IL23R";"IL2R";"IL4";"IL4R";"IL6";"IL6R";"IRAK";"JAK1";"JAK3";"NFAT";"RORgt";"SOCS1";"STAT1";"STAT3";"STAT4";"STAT5";"STAT6";"TCR";"TGFb";"TGFbR";"Tbet";];
# Global variables
points_to_be_steady = 10;
resolution = 0.0001;
time_slices = 100;
t_interval = 5;
t_maximal = 100;
# ODEs
function xdot = f(x,t)
    x1_2=1;
    x1_31=1;
    x1_35=1;
    x1_1=1;
    x1_26=1;
    x1_36=1;
    x1_decay=1;
    x2_2=1;
    x2_32=1;
    x2_36=1;
    x2_decay=1;
    x3_decay=1;
    x4_3=1;
    x4_decay=1;
    x5_25=1;
    x5_29=1;
    x5_30=1;
    x5_22=1;
    x5_36=1;
    x5_decay=1;
    x6_5=1;
    x6_decay=1;
    x7_2=1;
    x7_decay=1;
    x8_7=1;
    x8_decay=1;
    x9_decay=1;
    x10_32=1;
    x10_9=1;
    x10_decay=1;
    x11_26=1;
    x11_decay=1;
    x12_decay=1;
    x13_12=1;
    x13_32=1;
    x13_decay=1;
    x14_decay=1;
    x15_decay=1;
    x16_15=1;
    x16_decay=1;
    x17_14=1;
    x17_decay=1;
    x18_25=1;
    x18_2=1;
    x18_28=1;
    x18_decay=1;
    x19_27=1;
    x19_18=1;
    x19_decay=1;
    x20_2=1;
    x20_1=1;
    x20_36=1;
    x20_26=1;
    x20_decay=1;
    x21_20=1;
    x21_decay=1;
    x22_13=1;
    x22_decay=1;
    x23_27=1;
    x23_6=1;
    x23_decay=1;
    x24_21=1;
    x24_decay=1;
    x25_33=1;
    x25_decay=1;
    x26_2=1;
    x26_29=1;
    x26_24=1;
    x26_36=1;
    x26_35=1;
    x26_1=1;
    x26_26=1;
    x26_decay=1;
    x27_28=1;
    x27_36=1;
    x27_decay=1;
    x28_23=1;
    x28_4=1;
    x28_decay=1;
    x29_8=1;
    x29_16=1;
    x29_decay=1;
    x30_2=1;
    x30_10=1;
    x30_decay=1;
    x31_17=1;
    x31_decay=1;
    x32_19=1;
    x32_decay=1;
    x33_decay=1;
    x34_decay=1;
    x35_34=1;
    x35_decay=1;
    x36_2=1;
    x36_28=1;
    x36_36=1;
    x36_decay=1;
    xdot(1)=((-exp(0.5*50)+exp(-50*(((((1+x1_31+x1_35+x1_1)./(x1_31+x1_35+x1_1))*((x1_31*x(31)+x1_35*x(35)+x1_1*x(1))./(1+x1_31*x(31)+x1_35*x(35)+x1_1*x(1)))))*((1-((1+x1_2+x1_26+x1_36)./(x1_2+x1_26+x1_36))*((x1_2*x(2)+x1_26*x(26)+x1_36*x(36))./(1+x1_2*x(2)+x1_26*x(26)+x1_36*x(36)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x1_31+x1_35+x1_1)./(x1_31+x1_35+x1_1))*((x1_31*x(31)+x1_35*x(35)+x1_1*x(1))./(1+x1_31*x(31)+x1_35*x(35)+x1_1*x(1)))))*((1-((1+x1_2+x1_26+x1_36)./(x1_2+x1_26+x1_36))*((x1_2*x(2)+x1_26*x(26)+x1_36*x(36))./(1+x1_2*x(2)+x1_26*x(26)+x1_36*x(36)))))-0.5)))))-(x1_decay*x(1));
    xdot(2)=((-exp(0.5*50)+exp(-50*(((((1+x2_2+x2_32)./(x2_2+x2_32))*((x2_2*x(2)+x2_32*x(32))./(1+x2_2*x(2)+x2_32*x(32)))))*((1-((1+x2_36)./(x2_36))*((x2_36*x(36))./(1+x2_36*x(36)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x2_2+x2_32)./(x2_2+x2_32))*((x2_2*x(2)+x2_32*x(32))./(1+x2_2*x(2)+x2_32*x(32)))))*((1-((1+x2_36)./(x2_36))*((x2_36*x(36))./(1+x2_36*x(36)))))-0.5)))))-(x2_decay*x(2));
    xdot(3)=-(x3_decay*x(3));
    xdot(4)=((-exp(0.5*50)+exp(-50*((((1+x4_3)./(x4_3))*((x4_3*x(3))./(1+x4_3*x(3))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x4_3)./(x4_3))*((x4_3*x(3))./(1+x4_3*x(3))))-0.5)))))-(x4_decay*x(4));
    xdot(5)=((-exp(0.5*50)+exp(-50*(((((1+x5_25+x5_30+x5_22+x5_36)./(x5_25+x5_30+x5_22+x5_36))*((x5_25*x(25)+x5_30*x(30)+x5_22*x(22)+x5_36*x(36))./(1+x5_25*x(25)+x5_30*x(30)+x5_22*x(22)+x5_36*x(36)))))*((1-((1+x5_29)./(x5_29))*((x5_29*x(29))./(1+x5_29*x(29)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x5_25+x5_30+x5_22+x5_36)./(x5_25+x5_30+x5_22+x5_36))*((x5_25*x(25)+x5_30*x(30)+x5_22*x(22)+x5_36*x(36))./(1+x5_25*x(25)+x5_30*x(30)+x5_22*x(22)+x5_36*x(36)))))*((1-((1+x5_29)./(x5_29))*((x5_29*x(29))./(1+x5_29*x(29)))))-0.5)))))-(x5_decay*x(5));
    xdot(6)=((-exp(0.5*50)+exp(-50*((((1+x6_5)./(x6_5))*((x6_5*x(5))./(1+x6_5*x(5))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x6_5)./(x6_5))*((x6_5*x(5))./(1+x6_5*x(5))))-0.5)))))-(x6_decay*x(6));
    xdot(7)=((-exp(0.5*50)+exp(-50*((((1+x7_2)./(x7_2))*((x7_2*x(2))./(1+x7_2*x(2))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x7_2)./(x7_2))*((x7_2*x(2))./(1+x7_2*x(2))))-0.5)))))-(x7_decay*x(7));
    xdot(8)=((-exp(0.5*50)+exp(-50*((((1+x8_7)./(x8_7))*((x8_7*x(7))./(1+x8_7*x(7))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x8_7)./(x8_7))*((x8_7*x(7))./(1+x8_7*x(7))))-0.5)))))-(x8_decay*x(8));
    xdot(9)=-(x9_decay*x(9));
    xdot(10)=((-exp(0.5*50)+exp(-50*(((((1+x10_9)./(x10_9))*((x10_9*x(9))./(1+x10_9*x(9)))))*((1-((1+x10_32)./(x10_32))*((x10_32*x(32))./(1+x10_32*x(32)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x10_9)./(x10_9))*((x10_9*x(9))./(1+x10_9*x(9)))))*((1-((1+x10_32)./(x10_32))*((x10_32*x(32))./(1+x10_32*x(32)))))-0.5)))))-(x10_decay*x(10));
    xdot(11)=((-exp(0.5*50)+exp(-50*((((1+x11_26)./(x11_26))*((x11_26*x(26))./(1+x11_26*x(26))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x11_26)./(x11_26))*((x11_26*x(26))./(1+x11_26*x(26))))-0.5)))))-(x11_decay*x(11));
    xdot(12)=-(x12_decay*x(12));
    xdot(13)=((-exp(0.5*50)+exp(-50*(((((1+x13_12)./(x13_12))*((x13_12*x(12))./(1+x13_12*x(12)))))*((1-((1+x13_32)./(x13_32))*((x13_32*x(32))./(1+x13_32*x(32)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x13_12)./(x13_12))*((x13_12*x(12))./(1+x13_12*x(12)))))*((1-((1+x13_32)./(x13_32))*((x13_32*x(32))./(1+x13_32*x(32)))))-0.5)))))-(x13_decay*x(13));
    xdot(14)=-(x14_decay*x(14));
    xdot(15)=-(x15_decay*x(15));
    xdot(16)=((-exp(0.5*50)+exp(-50*((((1+x16_15)./(x16_15))*((x16_15*x(15))./(1+x16_15*x(15))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x16_15)./(x16_15))*((x16_15*x(15))./(1+x16_15*x(15))))-0.5)))))-(x16_decay*x(16));
    xdot(17)=((-exp(0.5*50)+exp(-50*((((1+x17_14)./(x17_14))*((x17_14*x(14))./(1+x17_14*x(14))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x17_14)./(x17_14))*((x17_14*x(14))./(1+x17_14*x(14))))-0.5)))))-(x17_decay*x(17));
    xdot(18)=((-exp(0.5*50)+exp(-50*(((((1+x18_25+x18_2)./(x18_25+x18_2))*((x18_25*x(25)+x18_2*x(2))./(1+x18_25*x(25)+x18_2*x(2)))))*((1-((1+x18_28)./(x18_28))*((x18_28*x(28))./(1+x18_28*x(28)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x18_25+x18_2)./(x18_25+x18_2))*((x18_25*x(25)+x18_2*x(2))./(1+x18_25*x(25)+x18_2*x(2)))))*((1-((1+x18_28)./(x18_28))*((x18_28*x(28))./(1+x18_28*x(28)))))-0.5)))))-(x18_decay*x(18));
    xdot(19)=((-exp(0.5*50)+exp(-50*(((((1+x19_18)./(x19_18))*((x19_18*x(18))./(1+x19_18*x(18)))))*((1-((1+x19_27)./(x19_27))*((x19_27*x(27))./(1+x19_27*x(27)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x19_18)./(x19_18))*((x19_18*x(18))./(1+x19_18*x(18)))))*((1-((1+x19_27)./(x19_27))*((x19_27*x(27))./(1+x19_27*x(27)))))-0.5)))))-(x19_decay*x(19));
    xdot(20)=((-exp(0.5*50)+exp(-50*(((((1+x20_26)./(x20_26))*((x20_26*x(26))./(1+x20_26*x(26)))))*((1-((1+x20_2+x20_1+x20_36)./(x20_2+x20_1+x20_36))*((x20_2*x(2)+x20_1*x(1)+x20_36*x(36))./(1+x20_2*x(2)+x20_1*x(1)+x20_36*x(36)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x20_26)./(x20_26))*((x20_26*x(26))./(1+x20_26*x(26)))))*((1-((1+x20_2+x20_1+x20_36)./(x20_2+x20_1+x20_36))*((x20_2*x(2)+x20_1*x(1)+x20_36*x(36))./(1+x20_2*x(2)+x20_1*x(1)+x20_36*x(36)))))-0.5)))))-(x20_decay*x(20));
    xdot(21)=((-exp(0.5*50)+exp(-50*((((1+x21_20)./(x21_20))*((x21_20*x(20))./(1+x21_20*x(20))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x21_20)./(x21_20))*((x21_20*x(20))./(1+x21_20*x(20))))-0.5)))))-(x21_decay*x(21));
    xdot(22)=((-exp(0.5*50)+exp(-50*((((1+x22_13)./(x22_13))*((x22_13*x(13))./(1+x22_13*x(13))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x22_13)./(x22_13))*((x22_13*x(13))./(1+x22_13*x(13))))-0.5)))))-(x22_decay*x(22));
    xdot(23)=((-exp(0.5*50)+exp(-50*(((((1+x23_6)./(x23_6))*((x23_6*x(6))./(1+x23_6*x(6)))))*((1-((1+x23_27)./(x23_27))*((x23_27*x(27))./(1+x23_27*x(27)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x23_6)./(x23_6))*((x23_6*x(6))./(1+x23_6*x(6)))))*((1-((1+x23_27)./(x23_27))*((x23_27*x(27))./(1+x23_27*x(27)))))-0.5)))))-(x23_decay*x(23));
    xdot(24)=((-exp(0.5*50)+exp(-50*((((1+x24_21)./(x24_21))*((x24_21*x(21))./(1+x24_21*x(21))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x24_21)./(x24_21))*((x24_21*x(21))./(1+x24_21*x(21))))-0.5)))))-(x24_decay*x(24));
    xdot(25)=((-exp(0.5*50)+exp(-50*((((1+x25_33)./(x25_33))*((x25_33*x(33))./(1+x25_33*x(33))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x25_33)./(x25_33))*((x25_33*x(33))./(1+x25_33*x(33))))-0.5)))))-(x25_decay*x(25));
    xdot(26)=((-exp(0.5*50)+exp(-50*(((((1+x26_29+x26_24+x26_35+x26_26)./(x26_29+x26_24+x26_35+x26_26))*((x26_29*x(29)+x26_24*x(24)+x26_35*x(35)+x26_26*x(26))./(1+x26_29*x(29)+x26_24*x(24)+x26_35*x(35)+x26_26*x(26)))))*((1-((1+x26_2+x26_36+x26_1)./(x26_2+x26_36+x26_1))*((x26_2*x(2)+x26_36*x(36)+x26_1*x(1))./(1+x26_2*x(2)+x26_36*x(36)+x26_1*x(1)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x26_29+x26_24+x26_35+x26_26)./(x26_29+x26_24+x26_35+x26_26))*((x26_29*x(29)+x26_24*x(24)+x26_35*x(35)+x26_26*x(26))./(1+x26_29*x(29)+x26_24*x(24)+x26_35*x(35)+x26_26*x(26)))))*((1-((1+x26_2+x26_36+x26_1)./(x26_2+x26_36+x26_1))*((x26_2*x(2)+x26_36*x(36)+x26_1*x(1))./(1+x26_2*x(2)+x26_36*x(36)+x26_1*x(1)))))-0.5)))))-(x26_decay*x(26));
    xdot(27)=((-exp(0.5*50)+exp(-50*((((1+x27_28+x27_36)./(x27_28+x27_36))*((x27_28*x(28)+x27_36*x(36))./(1+x27_28*x(28)+x27_36*x(36))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x27_28+x27_36)./(x27_28+x27_36))*((x27_28*x(28)+x27_36*x(36))./(1+x27_28*x(28)+x27_36*x(36))))-0.5)))))-(x27_decay*x(27));
    xdot(28)=((-exp(0.5*50)+exp(-50*((((1+x28_23+x28_4)./(x28_23+x28_4))*((x28_23*x(23)+x28_4*x(4))./(1+x28_23*x(23)+x28_4*x(4))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x28_23+x28_4)./(x28_23+x28_4))*((x28_23*x(23)+x28_4*x(4))./(1+x28_23*x(23)+x28_4*x(4))))-0.5)))))-(x28_decay*x(28));
    xdot(29)=((-exp(0.5*50)+exp(-50*((((1+x29_8+x29_16)./(x29_8+x29_16))*((x29_8*x(8)+x29_16*x(16))./(1+x29_8*x(8)+x29_16*x(16))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x29_8+x29_16)./(x29_8+x29_16))*((x29_8*x(8)+x29_16*x(16))./(1+x29_8*x(8)+x29_16*x(16))))-0.5)))))-(x29_decay*x(29));
    xdot(30)=((-exp(0.5*50)+exp(-50*(((((1+x30_10)./(x30_10))*((x30_10*x(10))./(1+x30_10*x(10)))))*((1-((1+x30_2)./(x30_2))*((x30_2*x(2))./(1+x30_2*x(2)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x30_10)./(x30_10))*((x30_10*x(10))./(1+x30_10*x(10)))))*((1-((1+x30_2)./(x30_2))*((x30_2*x(2))./(1+x30_2*x(2)))))-0.5)))))-(x30_decay*x(30));
    xdot(31)=((-exp(0.5*50)+exp(-50*((((1+x31_17)./(x31_17))*((x31_17*x(17))./(1+x31_17*x(17))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x31_17)./(x31_17))*((x31_17*x(17))./(1+x31_17*x(17))))-0.5)))))-(x31_decay*x(31));
    xdot(32)=((-exp(0.5*50)+exp(-50*((((1+x32_19)./(x32_19))*((x32_19*x(19))./(1+x32_19*x(19))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x32_19)./(x32_19))*((x32_19*x(19))./(1+x32_19*x(19))))-0.5)))))-(x32_decay*x(32));
    xdot(33)=-(x33_decay*x(33));
    xdot(34)=-(x34_decay*x(34));
    xdot(35)=((-exp(0.5*50)+exp(-50*((((1+x35_34)./(x35_34))*((x35_34*x(34))./(1+x35_34*x(34))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*((((1+x35_34)./(x35_34))*((x35_34*x(34))./(1+x35_34*x(34))))-0.5)))))-(x35_decay*x(35));
    xdot(36)=((-exp(0.5*50)+exp(-50*(((((1+x36_28+x36_36)./(x36_28+x36_36))*((x36_28*x(28)+x36_36*x(36))./(1+x36_28*x(28)+x36_36*x(36)))))*((1-((1+x36_2)./(x36_2))*((x36_2*x(2))./(1+x36_2*x(2)))))-0.5)))./((1-exp(0.5*50))*(1+exp(-50*(((((1+x36_28+x36_36)./(x36_28+x36_36))*((x36_28*x(28)+x36_36*x(36))./(1+x36_28*x(28)+x36_36*x(36)))))*((1-((1+x36_2)./(x36_2))*((x36_2*x(2))./(1+x36_2*x(2)))))-0.5)))))-(x36_decay*x(36));
endfunction
#Initial states:
#The following initial states correspond to all the steady states in the discrete system,
#you may use them (or any other you want) to evaluate the behaviour of the system.
#Remember that the nodes have the following order:
# "Foxp3";"GATA3";"IFNb";"IFNbR";"IFNg";"IFNgR";"IL10";"IL10R";"IL12";"IL12R";"IL17";"IL18";"IL18R";"IL2";"IL23";"IL23R";"IL2R";"IL4";"IL4R";"IL6";"IL6R";"IRAK";"JAK1";"JAK3";"NFAT";"RORgt";"SOCS1";"STAT1";"STAT3";"STAT4";"STAT5";"STAT6";"TCR";"TGFb";"TGFbR";"Tbet";
#initial_state = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
#initial_state = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
#initial_state = [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1];
#initial_state = [0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0];
#initial_state = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
initial_state = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
x0 = initial_state;
is_steady_state = false;
t_initial = 0;
t_final = t_initial + t_interval;
iteration_counter = 0;
while(is_steady_state==false & t_final<=t_maximal)
   t=linspace(t_initial,t_final,time_slices);
   x=lsode("f",x0,t);
   points_inside_resolution = 0;
   for i = (time_slices-points_to_be_steady+1 : time_slices)
        difference = abs(x(i,:) - x(i-1,:));
        if (all(all(difference<resolution)))
            points_inside_resolution++;
        endif
    endfor
    if(points_inside_resolution==points_to_be_steady)
        is_steady_state = true;
        steady_state = x(time_slices,:);
    endif
    final_time = t_final;
    x0 = x(time_slices,:);
    t_initial = t_final;
    t_final = t_initial + t_interval;
    iteration_counter++;
endwhile
   t=linspace(0,final_time,time_slices*iteration_counter);
   x=lsode("f",initial_state,t);

trans_t = t';
data = [trans_t x];
header = ["time Foxp3 GATA3 IFNb IFNbR IFNg IFNgR IL10 IL10R IL12 IL12R IL17 IL18 IL18R IL2 IL23 IL23R IL2R IL4 IL4R IL6 IL6R IRAK JAK1 JAK3 NFAT RORgt SOCS1 STAT1 STAT3 STAT4 STAT5 STAT6 TCR TGFb TGFbR Tbet"];
save -text - header data;