%%A function for generating the corresponding differential equations from a chemical reaction kinetic stoichiometric matrix
%%Created by wuqi
%%2023.12.19
%%revised at 2023.12.20
function dydt = odefunrevise(t,y)
reaction_num=3;
speice_num=3;
dydt = zeros(speice_num,1);
%y=zeros(speice_num);
%matrix_stnumber=zeros(speice_num,reaction_num);
%matrix_reactrate=zeros(speice_num,1);
matrix_stnumber_neg = -[1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,1,0;0,1,1,0,0,0,0,0,1,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,1,0,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,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,1,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0;0,0,0,1,1,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,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,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,1,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,0,0,0,0,0;0,0,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,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,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,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,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,1];%反应物化学计量数矩阵（前面加上负号）
matrix_stnumber_pos = [0,1,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,1;1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0;1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,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,0,0;0,0,0,2,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,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,0,0,0,0;0,0,0,1,1,1,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,0,0,0,0,0,0,0,0;0,0,0,0,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,1,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,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,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,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,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,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,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,1,0];  %生成物化学计量数矩阵
matrix_reactrate=[0.35 0.266e2 0.123e5 0.86e-3 0.82e-3 0.15e5 0.13e-3 0.24e5 0.165e5 0.9e4 0.22e-1 0.12e5 0.188e1 0.163e5 0.480e7 0.35e-3 0.175e-1 0.1e9 0.444e12 0.124e4 0.210e1 0.578e1 0.474e-1 0.178e4 0.312e1];%反应速率系数矩阵
sz = size(matrix_stnumber_pos);
matrix_source=zeros(sz);
for j=1:reaction_num
   index_neg = find(matrix_stnumber_neg(:,j)<0);
   index_pos = find(matrix_stnumber_pos(:,j)>0);
   product = 1
   for k=1:length(index_neg)
       product = product*y(index_neg(k))^abs(matrix_stnumber_neg(index_neg(k),j));
   end
   for k=1:length(index_neg)
       matrix_source(index_neg(k),j) = matrix_source(index_neg(k),j)+matrix_stnumber_neg(index_neg(k),j)*product;
   end
   for k=1:length(index_pos)
       matrix_source(index_pos(k),j) = matrix_source(index_pos(k),j)+matrix_stnumber_pos(index_pos(k),j)*product;
   end
   matrix_source(:,j) = matrix_reactrate(j)*matrix_source(:,j);
end
vector_source = sum(matrix_source, 2);
dydt = vector_source
end