求教MATLAB齿轮振动编程问题：求解八自由度二阶微分方程问题

冷风醉

如题：八自由度二阶微分方程问题。根据文献（附件1）编写的程序，不知道对不对？还望大神指点。谢谢。程序如下：

function dy=TorqueVibratory(t,y)
t
dy=zeros(8,2);
%齿轮轴承质量KG
m1=100; m2=100; g=9.81;
%转动惯量kg.m2
Jp=0.3; Jg=0.3; J1=1; J2=1;
%刚度系数
kt1=9000000; kt2=9000000; kx1=100000000; ky1=100000000; kx2=100000000; ky2=100000000;
%阻尼系数
ct1=400; ct2=400; cx1=500; cy1=500; cx2=500; cy2=500;
%主动轮转速
n1=1500; km=5*10^8;
omm1=n1*2*pi;
omm2=omm1;
%输入轴扭矩N.m
Tdm=400;  Tgm=Tdm; Tdr=300; Tgr=Tdr;
phid=omm1*t+y(1); phig=omm1*t+y(15);
Tp=Tdm+Tdr*sin(omm1*t+phid);
Tg=Tgm+Tgr*sin(omm2*t+phig);
%齿轮角度参数
alpha=20*pi/180;
%齿轮基础参数 单位mm
Z1=20; Z2=20; b1=16.5/1000; rb1=0.1;rb2=rb1; %n为齿轮转速
wm=2*pi*n1*Z1/60; em=2.0*10^-5; er=3.0*10^-5;  p1=1.5*10^-5; p2=p1; b=4.0*10^-5;%齿轮偏心距m
cm=1200;
%齿轮力计算

phi1=omm1*t+y(7); phi2=omm2*t+y(13);
et=em+er*sin(wm);%公式4
delta=(rb1*y(7)-rb2*y(13))+(y(3)-y(9)+p1*cos(-phi1)-p2*cos(phi2))*sin(alpha) ...
    +(y(5)-y(11)+p1*sin(-phi1)-p2*sin(phi2))*cos(alpha)-et;%公式3
deltat=(rb1*y(8)-rb2*y(14))+(y(4)-y(10)+p1*cos(-phi1)-p2*cos(phi2))*sin(alpha) ...
    +(y(6)-y(12)+p1*sin(-phi1)-p2*sin(phi2))*cos(alpha)-et;%公式3
if delta>b%公式6
    fdelta=delta-b;
else
    if delta<-b
        fdelta=delta+b;
    else
        fdelta=0;
    end
end
Fm=km*fdelta+cm*deltat;%公式5

%质量矩阵 p=ρ=齿轮偏心距;
M1=[Jp m1 m1 J1+m1*p1^2 m2 m2 J2+m2*p2^2 Jg];
M=diag(M1);

%刚度矩阵
K=zeros(8,8);
K(1,1)=kt1; K(1,4)=-kt1;
K(2,2)=kx1;
K(3,3)=ky1;
K(4,4)=ky1; K(4,1)=-kt1;
K(5,5)=kx2;
K(6,6)=ky2;
K(7,7)=kt2; K(7,8)=-kt2;
K(8,8)=kt2; K(8,7)=-kt2;

%阻尼矩阵 ω=ommiga=omm=转速
C=zeros(8,8);
C(1,1)=ct1; C(1,4)=-ct1;
C(2,2)=cx1;
C(3,3)=cy1;
C(4,4)=cy1; C(4,1)=-ct1;
C(5,5)=cx2;
C(6,6)=cy2;
C(7,7)=ct2; C(7,8)=-ct2;
C(8,8)=ct2; C(8,7)=-ct2;

%F矩阵 p=ρ=齿轮偏心距  Φ=phi  ω=ommiga=omm=转速 α=alpha
F=[Tp
    m1*p1*cos(omm1*t+y(7))*(omm1*t+y(8))^2+m1*p1*sin(omm1*t+y(7))*dy(8)-Fm*sin(alpha)
    -m1*g+m1*p1*sin(omm1*t+y(7))*(omm1*t+y(8))^2-m1*p1*cos(omm1*t+y(7))*dy(8)-Fm*cos(alpha)
    m1*p1*dy(4)*sin(omm1*t+y(7))-m1*p1*dy(6)*cos(omm1*t+y(7))-Fm*rb1
    m2*p2*cos(omm2*t+y(13))*(omm2*t+y(14))^2+m2*p2*sin(omm2*t+y(13))*dy(14)-Fm*sin(alpha)
    -m2*g+m2*p2*sin(omm2*t+y(13))*(omm2*t+y(14))^2-m2*p2*cos(omm2*t+y(13))*dy(14)-Fm*cos(alpha)
    m2*p2*dy(10)*sin(omm2*t+y(13))-m2*p2*dy(12)*cos(omm2*t+y(13))-Fm*rb2
    -Tg];
%系统动力学方程
dy=[y(2); (F(1,1)-C(1,1)*y(2)-C(1,4)*y(8)-K(1,1)*y(1)-K(1,4)*y(7))/M(1,1)
    y(4); (F(2,1)-C(2,2)*y(4)-K(2,2)*y(3))/M(2,2)
    y(6); (F(3,1)-C(3,3)*y(6)-K(3,3)*y(5))/M(3,3)
    y(8); (F(4,1)-C(4,4)*y(8)-C(4,1)*y(2)-K(4,4)*y(7)-K(4,1)*y(1))/M(4,4)
    y(10); (F(5,1)-C(5,5)*y(10)-K(5,5)*y(9))/M(5,5)
    y(12); (F(6,1)-C(6,6)*y(12)-K(6,6)*y(11))/M(6,6)
    y(14); (F(7,1)-C(7,7)*y(14)-C(7,8)*y(16)-K(7,7)*y(13)-K(7,8)*y(15))/M(7,7)
    y(16); (F(8,1)-C(8,8)*y(16)-C(8,7)*y(14)-K(8,8)*y(15)-K(8,8)*y(13))/M(8,8)];
end



调用程序如下：
clc;
clear;
t_span=linspace(4,4.05,20);
Y0=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
[t,y]=ode45('TorqueVibratory',t_span,Y0);
plot(t,y(:,4))
hold on

冷风醉

@dollfish000