t0=0
t=t0:0.05:100
s0=[1.0e-8;1.0e-8;1.0e-8;1.0e-8;1.0e-8;1.0e-8;1.0e-8;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;]
k1=5.0e6
k2=17
k3=3
k4=3
k6=1.0e3
k7=5.0e6
k8=1.0e6


function sfunc = sdot(t,s)
sfunc(1)=-k1*s(1)*s(2)-k1*s(1)*s(3)-k1*s(1)*s(22)-k1*s(1)*s(21)
sfunc(2)=-k1*s(1)*s(2)-k1*s(2)*s(9)-k1*s(2)*s(15)
sfunc(3)=-k1*s(1)*s(3)-k1*s(3)*s(8)-k1*s(3)*s(11)-k1*s(3)*s(13)-k1*s(3)*s(24)-k1*s(3)*s(28)-k1*s(3)*s(26)-k1*s(3)*s(31)-k1*s(3)*s(25)-k1*s(3)*s(35)
sfunc(4)=-k8*s(4)*s(15)
sfunc(5)=-k8*s(5)*s(35)

sfunc(6)=-k7*s(6)*s(15)
sfunc(7)=-k1*s(35)*s(7)
sfunc(8)=k1*s(1)*s(2)-k2*s(8)-k1*s(3)*s(8)-k1*s(8)*s(15)-k1*s(8)*s(21)
sfunc(9)=k1*s(1)*s(3)-k1*s(2)*s(9)-k1*s(9)*s(22)-k1*s(9)*s(21)
sfunc(10)=k1*s(3)*s(8)+k1*s(1)*s(9)-k2*s(10)-k1*s(10)*s(15)-k1*s(10)*s(21)

sfunc(11)=k2*s(8)-k1*s(3)*s(11)-k3*s(11)+k2*s(13)+k2*s(24)-k1*s(11)*s(21)
sfunc(12)=k2*s(10)+k1*s(3)*s(11)-k3*s(12)+k2*s(14)+k2*s(23)-k1*s(12)*s(21)
sfunc(13)=k3*s(11)-k1*s(3)*s(13)-k2*s(13)-k1*s(13)*s(21)
sfunc(14)=k3*s(12)+k1*s(3)*s(13)-k2*s(14)-k1*s(14)*s(21)
sfunc(15)=k2*s(13)+k2*s(14)-k7*s(6)*s(15)-k8*s(4)*s(15)-k1*s(8)*s(15)-k1*s(2)*s(15)+k2*s(23)-k1*s(10)*s(15)+k2*s(41)+k2*s(42)+k2*s(46)+k2*s(48)+k2*s(49)

sfunc(16)=k7*s(6)*s(15)-k2*s(16)
sfunc(17)=k2*s(16)
sfunc(18)=k2*s(16)
sfunc(19)=k8*s(4)*s(15)-k2*s(19)
sfunc(20)=k2*s(19)

sfunc(21)=k2*s(19)-k1*s(1)*s(21)-k1*s(11)*s(21)-k1*s(13)*s(21)-k1*s(9)*s(21)-k1*s(21)*s(23)-k1*s(21)*s(8)-k1*s(14)*s(21)-k1*s(10)*s(21)-k1*s(12)*s(21)-k1*s(24)*s(21)
sfunc(22)=k1*s(2)*s(15)-k1*s(1)*s(22)-k1*s(9)*s(22)
sfunc(23)=k1*s(3)*s(24)+k1*s(9)*s(22)-k2*s(23)+k1*s(10)*s(15)-k1*s(21)*s(23)
sfunc(24)=k1*s(8)*s(15)-k2*s(24)-k1*s(3)*s(24)+k1*s(1)*s(22)-k1*s(21)*s(24)
sfunc(25)=k1*s(1)*s(21)-k2*s(25)-k4*s(25)-k1*s(3)*s(25)

sfunc(26)=k4*s(31)+k4*s(25)-k6*s(26)-k1*s(3)*s(26)+k4*s(39)+k4*s(47)+k4*s(48)
sfunc(27)=k4*s(25)-k2*s(27)+k4*s(37)
sfunc(28)=-k1*s(3)*s(28)+k6*s(26)
sfunc(29)=k6*s(26)+k6*s(36)
sfunc(30)=k4*s(31)+k2*s(27)+k4*s(38)

sfunc(31)=k2*s(25)-k4*s(31)-k1*s(3)*s(31)
sfunc(32)=k1*s(3)*s(28)-k2*s(32)+k6*s(36)
sfunc(33)=k2*s(32)-k3*s(33)+k2*s(34)
sfunc(34)=k3*s(33)-k2*s(34)
sfunc(35)=k2*s(34)-k8*s(5)*s(35)-k1*s(35)*s(7)-k1*s(3)*s(35)

sfunc(36)=k1*s(3)*s(26)-k6*s(36)+k4*s(37)+k4*s(38)+k4*s(45)+k4*s(46)+k4*s(43)
sfunc(37)=k1*s(9)*s(21)+k1*s(3)*s(25)-k4*s(37)-k2*s(37)
sfunc(38)=k1*s(3)*s(31)+k2*s(37)-k4*s(38)
sfunc(39)=k1*s(11)*s(21)-k3*s(39)+k2*s(41)-k4*s(39)+k2*s(47)+k2*s(48)
sfunc(40)=k4*s(39)-k3*s(40)+k2*s(42)+k4*s(45)+k2*s(44)+k2*s(49)

sfunc(41)=k1*s(13)*s(21)+k3*s(39)-k2*s(41)
sfunc(42)=k3*s(40)-k2*s(42)+k4*s(46)
sfunc(43)=k1*s(21)*s(23)-k2*s(43)-k4*s(43)+k1*s(10)*s(21)
sfunc(44)=k4*s(43)-k2*s(44)+k4*s(47)
sfunc(45)=k2*s(43)-k4*s(45)-k3*s(45)+k2*s(46)+k1*s(12)*s(21)

sfunc(46)=k3*s(45)-k2*s(46)-k4*s(46)+k1*s(14)*s(21)
sfunc(47)=k1*s(8)*s(21)-k2*s(47)-k4*s(47)
sfunc(48)=k1*s(21)*s(24)-k2*s(48)-k4*s(48)
sfunc(49)=k4*s(48)-k2*s(49)
sfunc(50)=k8*s(5)*s(35)-k2*s(50)

sfunc(51)=k2*s(50)
sfunc(52)=k2*s(50)
sfunc(53)=k1*s(35)*s(7)-k2*s(53)
sfunc(54)=k2*s(53)
sfunc(55)=k2*s(53)


endfunction
s=ode("rkf",s0,t0,t,sdot)
plot2d(t,s(1,:),1)
plot2d(t,s(2,:),3)
plot2d(t,s(3,:),4)
//plot2d(t,s(4,:),1)
//plot2d(t,s(5,:),1)
//plot2d(t,s(6,:),1)
//plot2d(t,s(7,:),1)

//plot2d(t,s(15,:),3)
//plot2d(t,s(35,:),4)

plot2d(t,s(18,:),5)
plot2d(t,s(55,:),2)

temp=cat(1,t,s)

csvWrite(temp,"ABall.csv")