ode计算中可以使用元胞作为变量吗？

飞天踏雪

因为课题需要，我把变量考虑成元胞，程序如下
可调用不行，求高手指点
function myfuncircleflow1
c10=zeros(2,2);
c10(1,1)=1;
c10=c10(:);
c20=zeros(2,2);
c20=c20(:);
c30=zeros(2,2);
c30=c30(:);
c40=zeros(2,2);
c40=c40(:);
c50=zeros(2,2);
c50=c50(:);
c0={c10;c20;c30;c40;c50};
[t,y]=ode45(@myfuncircleflow,[0 8],c0);
plot(t,y);
function F=myfuncircleflow(t,c)
c1=size(4);
c2=size(4);
c3=size(4);
c4=size(4);
c5=size(4);

c={c1;c2;c3;c4;c5};


x=reshape(c1,2,2);
y=reshape(c2,2,2);
z=reshape(c3,2,2);
m=reshape(c4,2,2);
p=reshape(c5,2,2);



T=463;% K
k10=4.07e+6;%指前因子
k(1)=k10*exp(-7.88e+3/T);%速率常数
%i=1:2;
%j=1:2;

r1(1,1)=k(1)*x(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.0146);%反应速率
r1(2,1)=k(1)*x(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.0146);
r1(1,2)=k(1)*x(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.0146);
r1(2,2)=k(1)*x(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.0146);

k20=1.08e+6;
k(2)=k20*exp(-6.6e+3/T);
r2(1,1)=k(2)*y(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001)^0.5254;
r2(1,2)=k(2)*y(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001)^0.5254;
r2(2,1)=k(2)*y(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001)^0.5254;
r2(2,2)=k(2)*y(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001)^0.5254;


k30=9.8e+8;
k(3)=k30*exp(-11.16e+3/T);
r3(1,1)=k(3)*z(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001);
r3(1,2)=k(3)*z(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001);
r3(2,1)=k(3)*z(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001);
r3(2,2)=k(3)*z(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001);

k40=11.17e+9;
k(4)=k40*exp(-10.21e+3/T);
r4(1,1)=k(4)*m(1,1)/(1.427*x(1,1)+4.8419*m(1,1)+0.000001);
r4(1,2)=k(4)*m(1,2)/(1.427*x(1,2)+4.8419*m(1,2)+0.000001);
r4(2,1)=k(4)*m(2,1)/(1.427*x(2,1)+4.8419*m(2,1)+0.000001);
r4(2,2)=k(4)*m(2,2)/(1.427*x(2,2)+4.8419*m(2,2)+0.000001);

%----------------------------------------
%平衡方程
Q=0.0066;
q=0.0021;
%对PX
F=[Q*x(1,2)-Q*x(1,1)-r1(1,1);
    Q*x(1,1)-Q*x(2,1)-r1(2,1);
    Q*x(2,1)-Q*x(2,2)-r1(2,2);
    Q*x(2,2)-Q*x(1,2)-r1(1,2);
Q*y(1,2)-Q*y(1,1)+r1(1,1)-r2(1,1);
    Q*y(1,1)-Q*y(2,1)+r1(2,1)-r2(2,1);
    Q*y(2,1)-Q*y(2,2)+r1(2,2)-r2(2,2);
    Q*y(2,2)-Q*y(1,2)+r1(1,2)-r2(1,2);
Q*z(1,2)-Q*z(1,1)+r2(1,1)-r3(1,1);
    Q*z(1,1)-Q*z(2,1)+r2(2,1)-r3(2,1);
    Q*z(2,1)-Q*z(2,2)+r2(2,2)-r3(2,2);
    Q*z(2,2)-Q*z(1,2)+r2(1,2)-r3(1,2);
Q*m(1,2)-Q*m(1,1)+r3(1,1)-r4(1,1);
    Q*m(1,1)-Q*m(2,1)+r3(2,1)-r4(2,1);
    Q*m(2,1)-Q*m(2,2)+r3(2,2)-r4(2,2);
    Q*m(2,2)-Q*m(1,2)+r3(1,2)-r4(1,2);
Q*p(1,2)-Q*p(1,1)+r4(1,1);
    Q*p(1,1)-Q*p(2,1)+r4(2,1);
    Q*p(2,1)-Q*p(2,2)+r4(2,2);
    Q*p(2,2)-Q*p(1,2)+r4(1,2)];