t0=0
t=t0:0.1:6000
s0=[1.0e-7;1.0e-7;0;0;0.5e-7;0;0;1.0e-7;0;0;2e-7;0;0;0;0;0;0;0 ;0;0;0.5e-7;0;0;1.0e-7;0;0;2e-7;0;0;0;0;0;0;0 ;0;0;0.5e-7;0;0;1.0e-7;0;0;2e-7;0;0;0;0;0;0;0]

ks=5e4
kf=2e6
krf=10

function sfunc = sdot(t,s)
 sfunc(1) = -ks*s(1)*s(2)+ks*s(3)*s(4)
 sfunc(2) = -ks*s(1)*s(2)+ks*s(3)*s(4)
 sfunc(3) = ks*s(1)*s(2)-ks*s(3)*s(4)-kf*s(3)*s(11)+krf*s(18)-ks*s(3)*s(9)+ks*s(19)*s(20)
 sfunc(4) = ks*s(1)*s(2)-ks*s(3)*s(4)-kf*s(4)*s(5)-ks*s(4)*s(8)+ks*s(9)*s(10)+ks*s(10)*s(11)-ks*s(4)*s(12)-kf*s(4)*s(10)+krf*s(17)
 sfunc(5) = -kf*s(4)*s(5)
 sfunc(6) = kf*s(4)*s(5)
 sfunc(7) = kf*s(4)*s(5)
 sfunc(8) = -ks*s(4)*s(8)+ks*s(9)*s(10)-kf*s(8)*s(9)+krf*s(13)-kf*s(8)*s(11)+krf*s(16)
 sfunc(9) = ks*s(4)*s(8)-ks*s(9)*s(10)-kf*s(8)*s(9)+krf*s(13)-kf*s(9)*s(12)+krf*s(15)-ks*s(3)*s(9)+ks*s(19)*s(20)
 sfunc(10) = ks*s(4)*s(8)-ks*s(9)*s(10)-ks*s(10)*s(11)+ks*s(4)*s(12)-kf*s(4)*s(10)+krf*s(17)
 sfunc(11) = -ks*s(10)*s(11)+ks*s(4)*s(12)-kf*s(11)*s(12)+krf*s(14)-kf*s(8)*s(11)+krf*s(16)-kf*s(3)*s(11)+krf*s(18)
 sfunc(12) = -ks*s(4)*s(12)-ks*s(10)*s(11)-kf*s(11)*s(12)+krf*s(14)-kf*s(9)*s(12)+krf*s(15)
 sfunc(13) = kf*s(8)*s(9)-krf*s(13)
 sfunc(14) = kf*s(11)*s(12)-krf*s(14)
 sfunc(15) = kf*s(9)*s(12)-krf*s(15)
 sfunc(16) = kf*s(8)*s(11)-krf*s(16)
 sfunc(17) = kf*s(4)*s(10)-krf*s(17)
 sfunc(18) = kf*s(3)*s(11)-krf*s(18)

 sfunc(19) = ks*s(3)*s(9)-ks*s(19)*s(20)-kf*s(19)*s(27)+krf*s(34)-ks*s(19)*s(25)+ks*s(35)*s(36)
 sfunc(20) = ks*s(3)*s(9)-ks*s(19)*s(20)-kf*s(20)*s(21)-ks*s(20)*s(24)+ks*s(25)*s(26)+ks*s(26)*s(27)-ks*s(20)*s(28)-kf*s(20)*s(26)+krf*s(33)
 sfunc(21) = -kf*s(20)*s(21)
 sfunc(22) = kf*s(20)*s(21)
 sfunc(23) = kf*s(20)*s(21)
 sfunc(24) = -ks*s(20)*s(24)+ks*s(25)*s(26)-kf*s(24)*s(25)+krf*s(29)-kf*s(24)*s(27)+krf*s(32)
 sfunc(25) = ks*s(20)*s(24)-ks*s(25)*s(26)-kf*s(24)*s(25)+krf*s(29)-kf*s(25)*s(28)+krf*s(31)-ks*s(19)*s(25)+ks*s(35)*s(36)
 sfunc(26) = ks*s(20)*s(24)-ks*s(25)*s(26)-ks*s(26)*s(27)+ks*s(20)*s(28)-kf*s(20)*s(26)+krf*s(33)
 sfunc(27) = -ks*s(26)*s(27)+ks*s(20)*s(28)-kf*s(27)*s(28)+krf*s(30)-kf*s(24)*s(27)+krf*s(32)-kf*s(19)*s(27)+krf*s(34)
 sfunc(28) = -ks*s(20)*s(28)-ks*s(26)*s(27)-kf*s(27)*s(28)+krf*s(30)-kf*s(25)*s(28)+krf*s(31)
 sfunc(29) = kf*s(24)*s(25)-krf*s(29)
 sfunc(30) = kf*s(27)*s(28)-krf*s(30)
 sfunc(31) = kf*s(25)*s(28)-krf*s(31)
 sfunc(32) = kf*s(24)*s(27)-krf*s(32)
 sfunc(33) = kf*s(20)*s(26)-krf*s(33)
 sfunc(34) = kf*s(19)*s(27)-krf*s(34)

 sfunc(35) = ks*s(19)*s(25)-ks*s(35)*s(36)-kf*s(35)*s(43)+krf*s(50)
 sfunc(36) = ks*s(19)*s(25)-ks*s(35)*s(36)-kf*s(36)*s(37)-ks*s(36)*s(40)+ks*s(41)*s(42)+ks*s(42)*s(43)-ks*s(36)*s(44)-kf*s(36)*s(42)+krf*s(49)
 sfunc(37) = -kf*s(36)*s(37)
 sfunc(38) = kf*s(36)*s(37)
 sfunc(39) = kf*s(36)*s(37)
 sfunc(40) = -ks*s(36)*s(40)+ks*s(41)*s(42)-kf*s(40)*s(41)+krf*s(45)-kf*s(40)*s(43)+krf*s(48)
 sfunc(41) = ks*s(36)*s(40)-ks*s(41)*s(42)-kf*s(40)*s(41)+krf*s(45)-kf*s(41)*s(44)+krf*s(47)
 sfunc(42) = ks*s(36)*s(40)-ks*s(41)*s(42)-ks*s(42)*s(43)+ks*s(36)*s(44)-kf*s(36)*s(42)+krf*s(49)
 sfunc(43) = -ks*s(42)*s(43)+ks*s(36)*s(44)-kf*s(43)*s(44)+krf*s(46)-kf*s(40)*s(43)+krf*s(48)-kf*s(35)*s(43)+krf*s(50)
 sfunc(44) = -ks*s(36)*s(44)-ks*s(42)*s(43)-kf*s(43)*s(44)+krf*s(46)-kf*s(41)*s(44)+krf*s(47)
 sfunc(45) = kf*s(40)*s(41)-krf*s(45)
 sfunc(46) = kf*s(43)*s(44)-krf*s(46)
 sfunc(47) = kf*s(41)*s(44)-krf*s(47)
 sfunc(48) = kf*s(40)*s(43)-krf*s(48)
 sfunc(49) = kf*s(36)*s(42)-krf*s(49)
 sfunc(50) = kf*s(35)*s(43)-krf*s(50)

endfunction

s=ode("rkf",s0,t0,t,sdot)
plot2d(t,s(7,:),5)
plot2d(t,s(23,:),3)
plot2d(t,s(39,:),2)
plot2d(t,s(2,:),1)
plot2d(t,s(9,:),6)
plot2d(t,s(25,:),7)
plot2d(t,s(41,:),4)

temp=cat(1,t,s)
csvWrite(temp,"Enzyme-freeSystem1.csv")

stacksize('max')