[求助] 1stOpt编写常微分方程组计算反应速率问题已有2人参与

本人这么写代码：（[b]需要求每步反应速率k值（即为ka和kc），参数很多，不知能否计算出来；并且提供三组数据[/b]）
Parameters ka(10), kc1(10), kc2(10), kc3(10), kc4(10), kc5(10), kc6(10), kc7(10), kc8(10), kc9(10), kc10(10) ;       [b]这里的kc有方法简化下么？不然就要写10个了[/b]
Variable t, P(10), Q(10), F;   [b]P1~P10也可以这么定义变量么，直接用P(10)？[/b]
Plot t[x], P1, P2, Q1, Q2, F;   [b]不同组分的含量不一样，在一张图中会画出来看不出来？如果不同物质分开作图，怎么写？[/b]
Plot t[x], P3, P4, Q3, Q4;
ODEFunction
           P1'=-ka1*P1*F-P1*Sum(j=1:10)(kc1[j]*Q[j]);
           P2'=-ka2*P2*F-P2*Sum(j=1:10)(kc2[j]*Q[j])+Sum(j=1:(2-1))(kc2[j]*P2*Q[2-j]);   [b]式子中的数字2都和P2的2有关系，如果P2变成P3，那么其他的2也都为3，那么针对P3'到P10'，有没有更为方便的方法？下面再写一个P3'的，要一直写到P10'的
           P3'=-ka3*P3*F-P3*Sum(j=1:10)(kc3[j]*Q[j])+Sum(j=1:(3-1))(kc3[j]*P3*Q[3-j]);[/b]
           同样类似递推编写的方法也适用于Q1'到Q10'，因为都只是角标不一样，加和的物质产生变化而已，在matlab中就可以定义出来比较方便，针对Q1'如下：
           Q1'=-Q1*Sum(j=1:10)(kc1[j]*P[j])+ka1*P1*F;   [b]同理，Q2'=-Q2*Sum(j=1:10)(kc2[j]*P[j])+ka2*P2*F; 也要一直写到Q10'[/b]
           F'=-F*Sum(j=1:10)(ka[j]*P[j]);

Data t, P1, Q1, P2;
10, 30, 40, 50, 60, 75, 110, 150, 230;
10, 30, 40, 50, 60, 75, 110, 150, 230;
6.288097876, 4.91765691, 3.878352807, 3.298605506, 2.850889397, 2.241683203, 1.781706784, 0.910037534, 0.693987349;
0.257670552, 0.312407281, 0.244993608, 0.189709656, 0.133296065, 0.0803651, 0.033400217, 0.010654615, 0.005606342;
0.236493005, 0.607274759, 0.615323542, 0.623870854, 0.61234489, 0.562467845, 0.538393919, 0.347783178, 0.305460229;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 120, 170, 230, 305;
5.882896195, 4.593250328, 3.701714247, 2.910582239, 2.571710069, 2.416218035, 1.920578469, 1.68785903, 1.416030469, 1.287574203;
0.255794223, 0.297297319, 0.219253251, 0.161355927, 0.121348782, 0.081854014, 0.039094897, 0.016390052, 0.008676848, 0.005710866;
0.203988176, 0.534247028, 0.637609417, 0.634725031, 0.638037303, 0.569498973, 0.479824153, 0.477902906, 0.424218287, 0.411113563;
Data t, P1, Q1, P2;
10, 30, 50, 60, 70, 85, 130, 190, 237;
5.879605129, 4.718556629, 3.450891407, 3.045516702, 2.681262558, 2.428241655, 1.569209161, 1.234288027, 1.237674719;
0.22471112, 0.268468006, 0.177155403, 0.1332699, 0.094942111, 0.071844393, 0.019781466, 0.008520798, 0.00503865;
0.226510416, 0.569687824, 0.596106833, 0.599362637, 0.580304305, 0.666507432, 0.464066174, 0.389844442, 0.423282218;


最后再附上编写过的微分方程组部分的matlab程序，便于您对照。
function dydt=KineticsEqs(t, y, ka, kc)
m = 10;      
n = 10;
for i=1:m    %-----------------------将P赋值给y-----------------------------
    P(i) = y(i);
end
for j=1:n    %-----------------------将Q赋值给y-----------------------------
    Q(j) = y(m+j);
end
F = y(m+n+1);%-----------------------将F赋值给y-----------------------------
%----------------------------以下为动力学微分方程----------------------------
for j=1:n
    polynomial_1(j)=kc(1,j)*Q(j);
end
r(1) = -(  ka(1)*P(1)*F+P(1)*sum(polynomial_1)  ); %-------微分方程1--------

for i=1:m
    for j=1:n
    polynomial_2(j)=kc(i,j)*Q(j);
    end
    polysum_2(i)=sum(polynomial_2);
end
for i=2:m
    for j=1:i-1
         polynomial_3(j)=kc(i,j)*P(j)*Q(i-j);
    end
    polysum_3(i)=sum(polynomial_3);
end
r(i) = -(  ka(i)*P(i)*F+P(i)*polysum_2(i)-polysum_3(i)  );

for i=1:m
    for j=1:n
        polynomial_4(j)=kc(i,j)*P(j);
    end
    r(i+n) = -(  Q(i)*sum(polynomial_4)-ka(i)*P(i)*F  );%-----微分方程3-----
end

for j=1:n
        polynomial_5(j)=ka(j)*P(j);
end
r(m+n+1) = -(  F*sum(polynomial_5)  );             %-------微分方程4--------

dydt = r';