这个matlab程序运行没问题,但是不出结果，请高手指点

cxastudy

%  SFD—滑动轴承刚性转子系统的力学模型-位移响应
function Main

%转子系统各项参数
global  m g e omega c eta R L u E L1 I;  %笔记
m=10.2;              %转子质量
g=9.8;
e=0.0035;      %圆盘的偏心距
omega=300;          %omega为转子转速
c=0.5*10^(-5);      %转轴和轴承之间的半径间隙，c=轴承半径-转轴半径
eta=0.012;          %油膜的粘度系数
R=0.005;            %转轴半径
L=0.05;              %转轴长度
u=0.004;            %裂纹深度
E=2.1*10^11;        %弹性模量
L1=0.01;            %轴承宽度
I=pi*R^4/4;
t0=0:1:10;
x0=[0.0006 0.5 0.0008 0.7];
[t,zz]=ode45(@Weifen,t0,x0);         %显式常微分方程组的解法，第237页
n=length(t)-6;     %初始时间的长度减去150000，为了去掉瞬态响应
for i=1:6;                 %去掉前面的瞬态响应
    tt(i)=t(i+n);
    x(i)=zz(i+n,1);
    dx(i)=zz(i+n,2);     %位移的8个解都要去掉瞬态响应
    y(i)=zz(i+n,3);
    dy(i)=zz(i+n,4);
end
%输出图像，
figure(1)
plot(tt,x);
xlabel('时间/s');
ylabel('轴承中心x方向位移');
hold on   %图像保持，以下图片全可以用，第162页
figure(2)
plot(tt,y);
xlabel('时间/s');
ylabel('轴承中心y方向位移');
hold on

figure(3)
plot(x,y);
xlabel('轴承中心x方向位移');
ylabel('轴承中心y方向位移');
hold on

figure(4)
plot(x,dx);
xlabel('轴承中心x方向位移');
ylabel('轴承中心x方向速度');
hold on

function dz=Weifen(t0,z)                   %微分方程，function的用法见笔记
global  m e g omega;
%[epsilon1,depsilon1,phi1,dphi1]=transformation(z(1),z(3));
[Pe,Pt,sinj2,cosj2]=youmoli(z(1),z(2),z(3),z(4),omega) ;         %滑动轴承刚性转子系统轴承油膜力
[Kx,Ky,Kxy,Kyx]=huxigangdu(t0);
dz=zeros(4,1);
dz(1)=z(2);
dz(2)=e*omega^2*cos(omega*t0)+Pt*sinj2/m-Pe*cosj2/m-Kx*z(1)/m-Kxy*z(3)/m;
dz(3)=z(4);
dz(4)=e*omega^2*sin(omega*t0)-Pt*cosj2/m-Pe*sinj2/m-g-Ky*z(3)/m-Kyx*z(1)/m


function[Pe,Pt,sinj2,cosj2]=youmoli(x,dx,y,dy,omega)   %滑动轴承刚性转子系统轴承油膜力
%师兄pdf论文中56页公式
global c eta R L1;
h=sqrt(x^2+y ^2);
epsilon=h;
depsilon=x*y/(h^3);
%phi=asin(y/h)
%dphi=(y^2-x^2)/h^4
sinj2=y/h;
cosj2=x/h;
bb=x
vv=dx
dphi1=(dy*x-dx*y)/h^2;
F1=2*epsilon^2/((1-epsilon^2)*(2+epsilon^2));
F2=(pi/2-8/(pi*(2+epsilon^2)))/((1-epsilon^2)^(3/2));
F3=pi*epsilon/((1-epsilon^2)^(1/2)*(2+epsilon^2));
F4=2*epsilon/((1-epsilon^2)*(2+epsilon^2));
A=6*eta*R^3*L1/c^2;
Pe=A*((omega-2*dphi1)*F1+2*depsilon*F2)
Pt=A*((omega-2*dphi1)*F3+2*depsilon*F4)


function [Kx,Ky,Kxy,Kyx]=huxigangdu(t1)    %先不算裂纹位置
global u R E L omega I;
%syms t;
a=u/R;  %a的取值范围必须为0.5到1之间，否则运算会出现复数
I2x=pi*R^4/8-R^4/4*((1-a)*(1-4*a+2*a^2)*sqrt(2*a-a^2)+asin(1-a));
b=sqrt(a*(2-a));
Ar=R^2*(pi-acos(1-a)+(1-a)*b);
d=2/(3*Ar)*R^3*(2*a-a^2)^(3/2);
I11=I-I2x-Ar*d^2;
%gamma=R*sqrt(2*a-1);
I2y=R^4/12*((1-a)*(2*a^2-4*a-3)*b+3*asin(b));
%digits(10);
%I2y=vpa(I2y1);
I12=I-I2y;
syms k;
theta2=pi/2+acos(1-a);
theta12=0.8*theta2;
f=(2*theta12*sin(pi-k*theta12))/(pi^2-k^2*theta12^2);   %*sin(k*omega*t);
s=symsum(f,k,1,3);
digits(10);
ss=vpa(s);
Kx1=3*E/L^3*(I-(I-I11)*(1+cos(omega*t1)));
Ky1=3*E/L^3*(I-(I-I12)*(1+cos(omega*t1)));
Kxy1=3*E/L^3*(((I2x-I2y)/2+Ar*d^2/2)*ss*sin(omega*t1));
Kx=abs(vpa(Kx1));
Ky=abs(vpa(Ky1))
Kxy=abs(vpa(Kxy1))
Kyx=Kxy;