我列了一对齿的动力学方程，不知道哪错了，大家帮忙看看

姜洋

function dy=donglixue(t,y)
beta=18;              %螺旋角
p=560;                %输入功率，kw
n1=4100;              %输入转速，r/min
alfa=20;              %齿形角
d1=0.2081895;         %齿轮1分度圆直径，m
d2=0.5425545;         %齿轮2分度圆直径，m
g=0.08;               %啮合阻尼比
k=7850;               %铁的密度，kg/m3
km=18.553*10^9;       %平均刚度
i=2.6061;             %传动比
z1=33;                %齿轮1的齿数
b=0.065;              %齿宽，m
b3=2.5*10^(-4);       %侧隙
ea=1.08*10^(-4);      %误差幅值
f=pi/2;               %相位角
s=0.2;                %啮合刚度幅值ka/km
k1=0.5*10^6;          %齿轮1的y方向支承刚度
k2=0.5*10^6;          %齿轮1的z方向支承刚度
k3=0.5*10^6;          %齿轮2的y方向支承刚度
k4=0.5*10^6;          %齿轮2的z方向支承刚度
c1=0;                 %齿轮1的y方向支承阻尼
c2=0;                 %齿轮1的z方向支承阻尼
c3=0;                 %齿轮2的y方向支承阻尼
c4=0;                 %齿轮2的z方向支承阻尼
f1=0;                 %齿轮1支承系统所受的横向外载荷
f2=0;                 %齿轮2支承系统所受的横向外载荷
f3=0;                 %齿轮1支承系统所受的径向外载荷
f4=0;                 %齿轮2支承系统所受的径向外载荷
at=atan(tan(alfa*pi/180)/cos(beta*pi/180));%端面齿形角
r1=d1*cos(at)/2;      %齿轮1基圆半径，m
r2=d2*cos(at)/2;      %齿轮2基圆半径，m
m1=pi*(d1/2)^2*b*k;   %齿轮1的质量，kg
m2=pi*(d2/2)^2*b*k;   %齿轮2的质量，kg
j1=m1*r1^2/2;         %齿轮1的转动惯量
j2=m2*r2^2/2;         %齿轮2的转动惯量
me=j1*j2/(j1*r2^2+j2*r1^2);%质量当量
t1=9.55*10^6*p/n1;    %转矩，N.mm
n2=n1/i;              %输出转速，r/min
t2=9.55*10^6*p/n2;    %转矩，N.mm
wn=sqrt(km/me);       %固有频率
wh=(n1/60)*pi*z1/30;  %啮合频率
w=wh/wn;%激励频率
ch=2*g*sqrt(km*r1^2*r2^2*j1*j2/(r1^2*j1+r2^2*j2));%啮合阻尼
w1=sqrt(k1/m1);
w2=sqrt(k2/m1);
w3=sqrt(k3/m2);
w4=sqrt(k4/m2);
c11=c1/(2*m1*wn);
c13=ch*cos(beta*pi/180)/(2*m1*wn);
c22=c3/(2*m2*wn);
c23=ch*cos(beta*pi/180)/(2*m2*wn);
c34=c2/(2*m1*wn);
c36=ch*sin(beta*pi/180)/(2*m1*wn);
c45=c4/(2*m2*wn);
c46=ch*sin(beta*pi/180)/(2*m2*wn);
c53=c13;
c63=c23;
k11=w1^2/wn^2;
k13=(1+s*sin(w*t+f))*cos(beta*pi/180)*me/m1;
k22=w3^2/wn^2;
k23=(1+s*sin(w*t+f))*cos(beta*pi/180)*me/m2;
k34=w2^2/wn^2;
k36=(1+s*sin(w*t+f))*sin(beta*pi/180)*me/m1;
k45=w4^2/wn^2;
k46=(1+s*sin(w*t+f))*sin(beta*pi/180)*me/m2;
k53=k13;
k63=k23;
g1=f1/(m1*b3*wn^2);
g2=f2/(m2*b3*wn^2);
g3=f3/(m1*b3*wn^2);
g4=f4/(m2*b3*wn^2);
g5=t1*10^(-3)/(m1*r1*b3*wn^2);
g6=t2*10^(-3)/(m2*r2*b3*wn^2);
e=ea*sin(w*t+f)/b3;   %无量纲时变啮合误差
ey=e*cos(beta*pi/180);%y方向的无量纲时变啮合误差
ez=e*sin(beta*pi/180);%z方向的无量纲时变啮合误差
y1=y(1)+y(9)-y(3)+y(11)-ey;
y2=y(5)-y(7)+(-y(1)+y(2)-y(9)-y(11))*tan(beta*pi/180)-ez;
y11=y(2)+y(10)-y(4)+y(12)-w*ea*cos(beta*pi/180)*cos(w*t+f)/(b3*wn);
y22=y(6)-y(8)+(-y(2)+y(4)-y(10)-y(12))*tan(beta*pi/180)-w*ea*sin(beta*pi/180)*cos(w*t+f)/(b3*wn);
if y1>1
    fy1=y1-1;
else if abs(y1)<=1
        fy1=0;
    else
        fy1=y1+1;
    end
end
if y2>1
    fy2=y2-1;
else if abs(y2)<=1
        fy2=0;
    else
        fy2=y2+1;
    end
end
dy=[y(2);g1-2*c11*y(2)-k11*y(1)-k13*fy1-2*c13*y11;
    y(4);-g2-2*c22*y(4)-k22*y(3)+k23*fy1+2*c23*y11;
    y(6);g3-2*c34*y(6)-k34*y(5)-k36*fy2-2*c36*y22;
    y(8);-g4-2*c45*y(8)-k45*y(7)+k46*fy2+2*c46*y22;
    y(10);g5-k53*fy1-2*c53*y11;
    y(12);-g6+k63*fy1+2*c63*y11];

vib

程序表达的原公式最好也能贴出来，大家对比者比较容易看，直接看程序太费劲了

hhbhhy

最好多点说明。

xautmcf

建议先把你的想法和实现的意图说明一下!要不然大家看起来费劲!!