ode45本身是变步长的呀，怎么设置最小步长？我改用其它ode函数计算都会出现这个警告Warning: Failure at t=5.373174e-001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.908935e-015) at time t.，用ode4函数算了也出现这个问题，缩小tspan也不行，现贴上我的程序请高手指点。

clear

clc

miu=0.0178;

L=201.6e-3;

R=210e-3;

LD=L/(2*R);

C=0.699e-3;

pofai=C/R;

mg=104860;

m=mg/9.8;

ome=450;

pa=2*miu*ome/(pofai.^2);

Kxx =    6.453356740612708;

Kxy =   -4.133413328965316;

Kyx =1.146014052781378; Kyy = 7.846056548809618;

Dxx =1.145752328144106;Dxy =1.200350375824454;

Dyx =1.200350375824454;Dyy = 2.277131377326176;

KXxx = 4.310137875942596;KXxy =1.621686159880267;KXyy=-32.99946362265302;

DXxx = 84.83808962692011;DXxy =178.9053867925529;DXyy =122.8157676818844;

KYxx = 59.88546585149381;KYxy = 77.19598357923496;KYyy =18.20529335915087;

DYxx = 96.394530935465;DYxy = 237.0037140824726;DYyy = 222.0296129664884;

kxx=(pa*R*L/C)*Kxx*(C/(m*C*ome.^2));

kxy=(pa*R*L/C)*Kxy*(C/(m*C*ome.^2));

kyx=(pa*R*L/C)*Kyx*(C/(m*C*ome.^2));

kyy=(pa*R*L/C)*Kyy*(C/(m*C*ome.^2));

dxx=(pa*R*L/(C*ome))*Dxx*(C*ome/(m*C*ome.^2));

dxy=(pa*R*L/(C*ome))*Dxy*(C*ome/(m*C*ome.^2));

dyx=(pa*R*L/(C*ome))*Dyx*(C*ome/(m*C*ome.^2));

dyy=(pa*R*L/(C*ome))*Dyy*(C*ome/(m*C*ome.^2));

kxxx=(pa*R*L/(C*C))*KXxx*(C*C/(m*C*ome.^2));

kxxy=(pa*R*L/(C*C))*KXxy*(C*C/(m*C*ome.^2));

kxyy=(pa*R*L/(C*C))*KXyy*(C*C/(m*C*ome.^2));

dxxx=(pa*R*L/(C*C*ome))*DXxx*(C*C*ome/(m*C*ome.^2));

dxxy=(pa*R*L/(C*C*ome))*DXxy*(C*C*ome/(m*C*ome.^2));

dxyy=(pa*R*L/(C*C*ome))*DXyy*(C*C*ome/(m*C*ome.^2));

kyxx=(pa*R*L/(C*C))*KYxx*(C*C/(m*C*ome.^2));

kyxy=(pa*R*L/(C*C))*KYxy*(C*C/(m*C*ome.^2));

kyyy=(pa*R*L/(C*C))*KYyy*(C*C/(m*C*ome.^2));

dyxx=(pa*R*L/(C*C*ome))*DYxx*(C*C*ome/(m*C*ome.^2));

dyxy=(pa*R*L/(C*C*ome))*DYxy*(C*C*ome/(m*C*ome.^2));

dyyy=(pa*R*L/(C*C*ome))*DYyy*(C*C*ome/(m*C*ome.^2));

K=[kxx,kxy,kyx,kyy];

D=[dxx,dxy,dyx,dyy];

KX=[kxxx,dxxx,kxxy,dxxy,kxyy,dxyy];

KY=[kyxx,dyxx,kyxy,dyxy,kyyy,dyyy];

exy=0.02;

t_final=10000;

h_opt=odeset;

z0=zeros(4,1);

for i=1:4

z0(i)=1e-6;

end

[t,y]=ode45('Finite_Jeffcott_trance_unbalance',[0,t_final],z0,h_opt,K,D,KX,KY,exy);

function y=Finite_Jeffcott_trance_unbalance(t,z,flag,K,D,KX,KY,rou)

fx=D(1)*z(3)+D(2)*z(4)+K(1)*z(1)+K(2)*z(2)+KX(1)*z(1).^2+2*KX(2)*z(1)*z(3)+KX(3)*z(1)*z(2)+KX(4)*(z(3)*z(2)+z(1)*z(4))+KX(5)*z(2).^2+2*KX(6)*z(2)*z(4);

fy=D(3)*z(3)+D(4)*z(4)+K(3)*z(1)+K(4)*z(2)+KY(1)*z(1).^2+2*KY(2)*z(1)*z(3)+KY(3)*z(1)*z(2)+KY(4)*(z(3)*z(2)+z(1)*z(4))+KY(5)*z(2).^2+2*KY(6)*z(2)*z(4);

Tao=t；

y=[z(3);z(4);

-fx+rou*sin(Tao);

-fy+rou*cos(Tao)];

t

其中，z矩阵为所要求的未知数x和y组成的矩阵z=[x,y,x',y'];

急求高手指点，已经郁闷很久了。