参考程序：
p=[1 0.01 1 0;
   0 0.085 1 0.17;
   1 0.01 1 0;
   0 0.051 1 0.15;
   1 0.0067 70 0;
   1 0.15 0.21 0;
   0 1 130 0;
   1 0.01 0.1 0;
   1 0.01 0.0044 0];
WIJ=[1 0 1;
     2 1 1;
     2 9 -1;
     3 2 1;
     4 3 1;
     4 8 -1;
     5 4 1;
     6 5 1;
     6 7 -0.212;
     7 6 1;
     8 6 1;
     9 7 1];
n=9;
y0=1;
yt0=[0 0 0 0 0 0 0 0 0];
h=0.01;
T=0;
T0=0;
Tf=10;
nout=7;
A=diag(p(:,1));
B=diag(p(:,2));
C=diag(p(:,3));
D=diag(p(:,4));
m=length(WIJ(:,1));
w0=zeros(n,1);
w=zeros(n,n);
for k=1:m
     if (WIJ(k,2))==0;
         w0(WIJ(k,1))=WIJ(k,3);
     else w(WIJ(k,1),WIJ(k,2))=WIJ(k,3);
     end;
end;
Q=B-D*w;
Qn=inv(Q);
R=C*w-A;
V1=C*w0;
Ab=Qn*R;
b1=Qn*V1;
Y=yt0';
y=Y(nout);
t=T0;
N=round((Tf-T0)/h);
for i=1:N
     k1=Ab*Y+b1*y0;
     k2=Ab*(Y+h*k1/2)+b1*y0;
     k3=Ab*(Y+h*k2/2)+b1*y0;
     k4=Ab*(Y+h*k3)+b1*y0;
     Y=Y+h*(k1+2*k2+2*k3+k4)/6;
     y=[y,Y(nout)];
     T=[T,t+h];
     t=t+h;
end;
[T',y'];
plot(T,y);
参考结果：

用另一种思考方式：
nblocks=10;
n1=1;
d1=[0.01 1];
n2=[0.17 1];
d2=[0.085 0];
n3=1;
d3=[0.01 1];
n4=[0.15 1];
d4=[0.051 0];
n5=70;
d5=[0.0067 1];
n6=0.21;
d6=[0.15 1];
n7=130;
d7=[1 0];
n8=0.1;
d8=[0.01 1];
n9=0.0044;
d9=[0.01 1];
n10=0.212;
d10=1;
blkbuild;
sys=ss(a,b,c,d);
q=[1 0 0;
   2 1 -9;
   3 2 0;
   4 3 -8;
   5 4 0;
   6 5 -10;
   7 6 0;
   8 6 0;
   9 7 0;
   10 7 0];
inputs=1;
outputs=7;
sysc=connect(sys,q,inputs,outputs)
step(tf(sysc))
参考结果：

两种方法结果比较(法一用红线表示，法二用蓝线表示)：

Simulink方法：

────────────
　　北京信息科技大学自动化学院“控制系统数字仿真与CAD”实验内容，仅供参考

08.05.26最后更新