一个程序无法调通，很急，哪位兄弟姐妹能帮我？我出车费和酬谢

sniper72

这是我的一个程序，一个制冷机的动态模型，没法调通，哪位兄弟姐妹可以帮助我的，我愿意出车费并有酬谢，我住上海杨浦区。很急，希望哪位能救救急。
电话：13818088246
kg=1;
fg=0.68;
kk=4;
fk=0.67;
ke=2.3;
fe=0.62;
ka=1;
fa=1.46;
kex=1;
fex=0.17;
d=0.0015;
te=10;
t2=35;
t4=80;
t7=65;
t8=50;
t9=40;
twmid=35;
twout=38;
tcout=13;
thout=92;
tk=40;
T=[te t2 t4 t7 t8 t9 twmid twout tcout thout tk]';
for n=1:500
    if n<200
        mst0=100;
        msw0=600;
        mw0=1.5;
        mh0=1.3;
        mc0=0.7;
        thin=95;
        tcin=18;
        twin=32;
    else
        mst0=160;
        msw0=600;
        mw0=1.5;
        mh0=1.3;
        mc0=0.7;
        thin=95;
        tcin=18;
        twin=32;
    end
mst=mst0*1.7/3600;
msw=msw0*1.7/3600;
mw=mw0/3.6;
mh=mh0/3.6;
mc=mc0/3.6;
cps=2;
cpw=4.2;
cm=0.46;  %钢的热容
%钢的质量
mg=50;
ma=100;
mk=40;
me=50;
mex=30;
re=2475;
rk=2427;
%工质在各单元中的停留时间
timeg=15;
timea=60;
timek=15;
timee=60;
timeex=15;
%外界热源工质在单元中停留时间
timeh=5;
timewk=5;
timewa=5;
timec=5;
U=[thin tcin twin d]';
A=[0 0 -mst*cps-0.65*kg*fg mst*cps-0.35*kg*fg 0 0 0 0 0 0.35*kg*fg 0;...
        0 0 0.65*kg*fg 0.35*kg*fg 0 0 0 0 0 -mh*cpw-0.5*kg*fg 0;...
        0 0 0 0 0 0 0.35*kk*fk 0.65*kk*fk 0 0 -kk*fk;...
        0 0 0 0 0 0 mw*cpw-0.35*kk*fk -mw*cpw-0.65*kk*fk 0 0 kk*fk;...
        0 -msw*cps-0.65*ka*fa 0 0 0 msw*cps-0.35*ka*fa -0.35*ka*fa 0 0 0 0;...
        0 0.65*ka*fa 0 0 0 0.35*ka*fa -mw*cpw-0.35*ka*fa 0 0 0 0;...
        -ke*fe 0 0 0 0 0 0 0 0.35*ke*fe 0 0;...
        ke*fe 0 0 0 0 0 0 0 -mc*cpw-0.35*ke*fe 0 0;...
        0 mst*cps-0.65*kex*fex 0.35*kex*fex -mst*cps-0.35*kex*fex 0.65*kex*fex 0 0 0 0 0 0;...
        0 0.65*kex*fex mst*cps-0.35*kex*fex 0.35*kex*fex -mst*cps-0.65*kex*fex 0 0 0 0 0 0;...
        0 (msw-mst)*cps 0 0 0 0 0 mst*cps 0 0 0];
B=[0.65*kg*fg 0 0 -rk;...
        mh*cpw-0.5*kg*fg 0 0 0;...
        0 0 0 rk;...
        0 0 0 0;...
        0 0 0.65*ka*fa re;...
        0 0 mw*cpw-0.65*ka*fa 0;...
        0 0.65*ke*fe 0 -re;...
        0 mc*cpw-0.65*ke*fe 0 0;...
        0 0 0 0;...
        0 0 0 0;...
        0 0 0 0];
Z=[0 0 1 1 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 1 1 0 0 0;...
        0 1 0 0 0 1 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 0 0 0 1 0 0;...
        0 1 0 1 0 0 0 0 0 0 0;...
        0 0 1 0 1 0 0 0 0 0 0;...
        0 0 0 0 0 1 0 0 0 0 0];
J=[(mst*timeg*cps+mg*cm)/2 0 0 0 0 0 0 0 0 0 0;...
        0 mh*timeh*cpw 0 0 0 0 0 0 0 0 0;...
        0 0 mk*cm 0 0 0 0 0 0 0 0;...
        0 0 0 mw*timewk*cpw 0 0 0 0 0 0 0;...
        0 0 0 0 (ma*cm+msw*timea*cps)/2 0 0 0 0 0 0;...
        0 0 0 0 0 mw*timewa*cpw 0 0 0 0 0;...
        0 0 0 0 0 0 me*cm+20*cpw 0 0 0 0;...
        0 0 0 0 0 0 0 mc*timec*cpw 0 0 0;...
        0 0 0 0 0 0 0 0 (mst*timeex*cps+mex*cm)/2 0 0;...
        0 0 0 0 0 0 0 0 0 (mst*timeex*cps+mex*cm)/2 0;...
        0 0 0 0 0 0 0 0 0 0 msw];
% T=[te t2 t4 t7 t8 t9 twmid twout tcout thout tk]
C=Z\(J\A);
D=Z\(J\B);
% G=expm(C)
%T=(C+eye)*T+D*U;
T=(C+eye(11,11))*T+D*U;
te=T(1);
t2=T(2);
t4=T(3);
t7=T(4);
t8=T(5);
t9=T(6);
twmid=T(7);
twout=T(8);
tcout=T(9);
thout=T(10);
tk=T(11);
d=(0.3614+0.4974*t4-0.6224*tk+1.2109e-4*t4*tk-0.5448*t2+1.7715*te-0.0299*t2*te)*mst/100;
end