matlab ode45求解齿轮动力学,使用ODE45求解齿轮系统动力学方程后结果发散

function [dx,ff1,ff2]=myfun(t,x)

t

beita=26;

mn=0.004;

z1=46;

z2=43;

z3=122;

T_in=200;

T_out=80;

roug1=7.8E3;

roug2=7.8E3;

roug3=7.8E3;

alphan=20;

alphat=atand(tand(alphan)/cosd(beita));

d1=z1*mn/cosd(beita)/1000;

db1=d1*cosd(alphat)/1000;

d2=z2*mn/cosd(beita)/1000;

db2=d2*cosd(alphat)/1000;

d3=z3*mn/cosd(beita)/1000;

db3=d3*cosd(alphat)/1000;

bp1=116/1000;

bp2=116/1000;

bp3=116/1000;

bp=116/1000;

I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;

I2=((roug2*pi*(d2/2)^2*bp2)*(d2/2)^2)/2;

I3=((roug3*pi*(d3/2)^2*bp3)*(d1/2)^2)/2;

m1=roug1*pi*(d1/2)^2*bp1;

m2=roug2*pi*(d2/2)^2*bp2;

m3=roug3*pi*((d3)-(d1+d2))^2*bp3;

r1=d1*cosd(alphat)/2;

r2=d2*cosd(alphat)/2;

r3=d3*cosd(alphat)/2;

fai_sp1x=90;

fai_s