function dy=order22(t,y);%用积分形式定义函数；

dy=zeros(16,1);

mc=854.97;mp=745.26;ms=332.69;m1=2135.1;m2=94.23;m3=580.72;m4=37.68;%质量赋值；

a=22.5*pi/180;ksp1=4.3313E9;ksp2=4.313e9;ksp3=4.3313E9;krp1=3.2989E9;krp2=3.2989E9;krp3=3.2989E9;ks1=5.045E7;K34=2.4276E9;K12=2.7783e9;kc=6E11;%刚度赋值

rbs=174.782e-3;rb1=503.2e-3;rb2=105.71e-3;rb3=348.625e-3;rbc=704.69e-3;rb4=88.801e-3;%基圆半径赋值；

cs1=2.476E3;csp=6.3211E3;crp=25.657E3;C12=1.6547E3;C34=0.3024E3;%阻尼赋值；

Tin=8.22E6+8.22E6*sin((38*pi)*t);Tout=Tin/94;%输入输出转矩；

F=[Tin/rbc,0,0,0,0,0,0,Tout/rb4];%F矩阵；

M_1=[mc+mp+mp+mp,mp,mp,mp,ms,m1,m2+m3*(rb3^2)/(rb2^2),m4];%质量矩阵

M=diag(M_1);

k11=[kc+(krp1+ksp1)*(cos(a)^2)+(krp2+ksp2)*(cos(a)^2)+(krp3+ksp3)*(cos(a)^2)];

k12=[(krp1-ksp1)*cos(a)];

k21=[(krp1-ksp1)*cos(a)];

k31=[(krp2-ksp2)*cos(a)];

k41=[(krp3-ksp3)*cos(a)];

k13=[(krp2-ksp2)*cos(a)];

k14=[(krp3-ksp3)*cos(a)];

k22=krp1+ksp1;

k33=krp2+ksp2;

k15=[-ksp1*cos(a)-ksp2*cos(a)-ksp3*cos(a)];

k51=k15;

k44=krp3+ksp3;

k55=[ksp1+ksp2+ksp3+ks1/(rbs)^2];

k56=[-ks1/(rbs*rb1)];

k65=k56;

k66=[K12+ks1/(rb1)^2];

k77=[K12+K34*(rb3)^2/(rb2)^2];

k78=[K34*rb3/rb2];

k87=k78;

k=[k11,k12,k13,k14,k15,0,0,0;

k21,k22,0,0,ksp1,0,0,0;

k31,0,k33,0,ksp2,0,0,0;

k41,0,0,k44,ksp3,0,0,0;

k51,ksp1,ksp2,ksp3,k55,k56,0,0;

0,0,0,0,k65,k66,K12,0;

0,0,0,0,0,K12,k77,k78;

0,0,0,0,0,0,k87,K34];%刚度矩阵；

c11=3*(crp+csp)*(cos(a)^2);

c12=[(crp-csp)*cos(a)];

c13=[(crp-csp)*cos(a)];

c14=[(crp-csp)*cos(a)];

c15=-3*csp*cos(a);

c21=[(crp-csp)*cos(a)];

c22=crp+csp;

c25=csp;

c31=[(crp-csp)*cos(a)];

c33=crp+csp;

c35=csp;

c41=c14;

c44=c33;

c45=c35;

c51=c15;

c52=c35;

c53=c35;

c54=c35;

c55=3*csp+cs1/(rbs)^2;

c56=-cs1/(rbs*rb1);

c65=c56;

c66=C12+cs1/(rb1)^2;

c67=C12;

c76=c67;

c77=C12+C34*(rb3)^2/(rb2)^2;

c78=C34*rb3/rb2;

c87=c78;

c88=C34;

c=[c11,c12,c13,c14,c15,0,0,0;

c21,c22,0,0,c25,0,0,0;

c31,0,c33,0,c35,0,0,0;

c41,0,0,c44,c45,0,0,0;

c51,c52,c53,c54,c55,c56,0,0;

0,0,0,0,c65,c66,c67,0;

0,0,0,0,0,c76,c77,c78;

0,0,0,0,0,0,c87,c88];%阻尼矩阵

M1=inv(M);

[cx, cy] = size(c);

[M1x, M1y] = size(M1);

for i=1:cx

for j=1:M1y

w11(i,j) =-M1(i,j)*c(i,j);

end                        %  -(M的逆矩阵)*C矩阵；

end

[kx, ky] = size(k);

[M1x, M1y] = size(M1);

for i=1:kx

for j=1:M1y

w12(i,j) =-M1(i,j)*k(i,j);

end

end                       %  -(M的逆矩阵)*K矩阵；

q1=eye(8);

q2=zeros(8,8);

q3=zeros(8,1);

V1=zeros(8,1);

V1(1)=M(1,1)*F(1)+M(1,2)*F(2)+M(1,3)*F(3)+M(1,4)*F(4)+M(1,5)*F(5)+M(1,6)*F(6)+M(1,7)*F(7)+M(1,8)*F(8);

V1(2)=M(2,1)*F(1)+M(2,2)*F(2)+M(2,3)*F(3)+M(2,4)*F(4)+M(2,5)*F(5)+M(2,6)*F(6)+M(2,7)*F(7)+M(2,8)*F(8);

V1(3)=M(3,1)*F(1)+M(3,2)*F(2)+M(3,3)*F(3)+M(3,4)*F(4)+M(3,5)*F(5)+M(3,6)*F(6)+M(3,7)*F(7)+M(3,8)*F(8);

V1(4)=M(4,1)*F(1)+M(4,2)*F(2)+M(4,3)*F(3)+M(4,4)*F(4)+M(4,5)*F(5)+M(4,6)*F(6)+M(4,7)*F(7)+M(4,8)*F(8);

V1(5)=M(5,1)*F(1)+M(5,2)*F(2)+M(5,3)*F(3)+M(5,4)*F(4)+M(5,5)*F(5)+M(5,6)*F(6)+M(5,7)*F(7)+M(5,8)*F(8);

V1(6)=M(6,1)*F(1)+M(6,2)*F(2)+M(6,3)*F(3)+M(6,4)*F(4)+M(6,5)*F(5)+M(6,6)*F(6)+M(6,7)*F(7)+M(6,8)*F(8);

V1(7)=M(7,1)*F(1)+M(7,2)*F(2)+M(7,3)*F(3)+M(7,4)*F(4)+M(7,5)*F(5)+M(7,6)*F(6)+M(7,7)*F(7)+M(7,8)*F(8);

V1(8)=M(8,1)*F(1)+M(8,2)*F(2)+M(8,3)*F(3)+M(8,4)*F(4)+M(8,5)*F(5)+M(8,6)*F(6)+M(8,7)*F(7)+M(8,8)*F(8);   %  M矩阵  *  K矩阵(不太会使用循环，说以逐一写的语句)

W1=[w11,w12];

W2=[q1,q2];

W=[W1;W2];

V=[V1;q3];%矩阵合并

dy(1)=W(1,1)*y(1)+W(1,2)*y(2)+W(1,3)*y(3)+W(1,4)*y(4)+W(1,5)*y(5)+W(1,6)*y(6)+W(1,7)*y(7)+W(1,8)*y(8)+V(1);

dy(2)=W(2,1)*y(1)+W(2,2)*y(2)+W(2,3)*y(3)+W(2,4)*y(4)+W(2,5)*y(5)+W(2,6)*y(6)+W(2,7)*y(7)+W(2,8)*y(8)+V(2);

dy(3)=W(3,1)*y(1)+W(3,2)*y(2)+W(3,3)*y(3)+W(3,4)*y(4)+W(3,5)*y(5)+W(3,6)*y(6)+W(3,7)*y(7)+W(3,8)*y(8)+V(3);

dy(4)=W(4,1)*y(1)+W(4,2)*y(2)+W(4,3)*y(3)+W(4,4)*y(4)+W(4,5)*y(5)+W(4,6)*y(6)+W(4,7)*y(7)+W(4,8)*y(8)+V(4);

dy(5)=W(5,1)*y(1)+W(5,2)*y(2)+W(5,3)*y(3)+W(5,4)*y(4)+W(5,5)*y(5)+W(5,6)*y(6)+W(5,7)*y(7)+W(5,8)*y(8)+V(5);

dy(6)=W(6,1)*y(1)+W(6,2)*y(2)+W(6,3)*y(3)+W(6,4)*y(4)+W(6,5)*y(5)+W(6,6)*y(6)+W(6,7)*y(7)+W(6,8)*y(8)+V(6);

dy(7)=W(7,1)*y(1)+W(7,2)*y(2)+W(7,3)*y(3)+W(7,4)*y(4)+W(7,5)*y(5)+W(7,6)*y(6)+W(7,7)*y(7)+W(7,8)*y(8)+V(7);

dy(8)=W(8,1)*y(1)+W(8,2)*y(2)+W(8,3)*y(3)+W(8,4)*y(4)+W(8,5)*y(5)+W(8,6)*y(6)+W(8,7)*y(7)+W(8,8)*y(8)+V(8);

dy(9)=W(9,1)*y(1)+W(9,2)*y(2)+W(9,3)*y(3)+W(9,4)*y(4)+W(9,5)*y(5)+W(9,6)*y(6)+W(9,7)*y(7)+W(9,8)*y(8)+V(9);

dy(10)=W(10,1)*y(1)+W(10,2)*y(2)+W(10,3)*y(3)+W(10,4)*y(4)+W(10,5)*y(5)+W(10,6)*y(6)+W(10,7)*y(7)+W(10,8)*y(8)+V(10);

dy(11)=W(11,1)*y(1)+W(11,2)*y(2)+W(11,3)*y(3)+W(11,4)*y(4)+W(11,5)*y(5)+W(11,6)*y(6)+W(11,7)*y(7)+W(11,8)*y(8)+V(11 );

dy(12)=W(12,1)*y(1)+W(12,2)*y(2)+W(12,3)*y(3)+W(12,4)*y(4)+W(12,5)*y(5)+W(12,6)*y(6)+W(12,7)*y(7)+W(12,8)*y(8)+V(12);

dy(13)=W(13,1)*y(1)+W(13,2)*y(2)+W(13,3)*y(3)+W(13,4)*y(4)+W(13,5)*y(5)+W(13,6)*y(6)+W(13,7)*y(7)+W(13,8)*y(8)+V(13);

dy(14)=W(14,1)*y(1)+W(14,2)*y(2)+W(14,3)*y(3)+W(14,4)*y(4)+W(14,5)*y(5)+W(14,6)*y(6)+W(14,7)*y(7)+W(14,8)*y(8)+V(14);

dy(15)=W(15,1)*y(1)+W(15,2)*y(2)+W(15,3)*y(3)+W(15,4)*y(4)+W(15,5)*y(5)+W(15,6)*y(6)+W(15,7)*y(7)+W(15,8)*y(8)+V(15);

dy(16)=W(16,1)*y(1)+W(16,2)*y(2)+W(16,3)*y(3)+W(16,4)*y(4)+W(16,5)*y(5)+W(16,6)*y(6)+W(16,7)*y(7)+W(16,8)*y(8)+V(16);% dy=W*y+V;

end

y0=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];tspan=[0 10];[t,y]=ode45(@order22,tspan,y0);

我的程序和你的类似，有的地方比较繁琐，跪求指点啊，困扰很久了