希望各位耐心赐教，谢谢。
% 函数文件
function dy=rigid(t,y)

global w
% 参数取值
RT=0.5;
RB=0.4;
sita1=340*pi/180;
sita2=290*pi/180;
sita3=250*pi/180;
sita4=200*pi/180;
sita5=160*pi/180;
sita6=110*pi/180;
sita7=70*pi/180;
sita8=20*pi/180;
A=(RT^2-RB^2)/2;
T1=cos(sita1)-cos(sita2);
T2=cos(sita3)-cos(sita4);
T3=cos(sita5)-cos(sita6);
T4=cos(sita7)-cos(sita8);
E1=sin(sita1)-sin(sita2);
E2=sin(sita3)-sin(sita4);
E3=sin(sita5)-sin(sita6);
E4=sin(sita7)-sin(sita8);
V=300;
P=50;
b1=25.42*pi/180;
b2=20.42*pi/180;
K1=0.82;
C=V^2*sin(b1)*P*(cos(b1)+K1*cos(b2));

h=0.5;
FX=A*C*T1+A*C*T2+A*C*T3+h*A*C*T4;
FY=A*C*E1+A*C*E2+A*C*E3+h*A*C*E4;

dy=zeros(8,1);
u=100;
cz=0.0002;
L=0.0285;
r=0.057;
B=atan((y(3)+2*y(4))/(y(1)-2*y(4)))-pi/2*sign((y(3)+2*y(4))/(y(1)-2*y(4)))-pi/2*sign(y(3)+2*y(2));
G=2*(pi/2+atan((y(3)*cos(B)-y(1)*sin(B))/(1-y(1)^2-y(3)^2)^0.5))/(1-y(1)^2-y(3)^2)^0.5;
V=(2+(y(3)*cos(B))-y(1)*sin(B)*G)/(1-(y(1)^2+y(3)^2));
S=(y(1)*cos(B)+y(3)*sin(B))/(1-(y(1)*cos(B)+y(3)*sin(B))^2);
o=u*w*r*L*(r/cz)^2*(L/(2*r))^2;
fx=-o*((y(1)-2*y(4))^2+(y(3)+2*y(2))^2)^0.5*(3*y(1)*V-sin(B)*G-2*cos(B)*S)/(1-y(1)^2-y(3)^2);
fy=-o*((y(1)-2*y(4))^2+(y(3)+2*y(2))^2)^0.5*(3*y(3)*V-cos(B)*G-2*sin(B)*S)/(1-y(1)^2-y(3)^2);

c1=9300;
c2=9300;
m1=50;
m2=420;
k=1.052*10^8;
g=9.8;
e=0.0004;
% 微分方程组
dy(1)=y(2);
dy(2)=-c1*y(2)/(m1*w)-k/2*(y(1)-y(5))/(m1*w^2)+fx/(m1*cz*w^2);

dy(3)=y(4);
dy(4)=-c1*y(4)/(m1*w)-k/2*(y(3)-y(7))/(m1*w^2)+fy/(m1*cz*w^2)-g/(w^2*cz);

dy(5)=y(6);
dy(6)=-c2*y(6)/(m2*w)-k*(y(5)-y(1))/(m2*w^2)+FX/(m2*cz*w^2)+e*cos(t)/cz;

dy(7)=y(8);
dy(8)=-c2*y(8)/(m2*w)-k*(y(7)-y(3))/(m2*w^2)+FY/(m2*cz*w^2)+e*sin(t)/cz-g/(cz*w^2);


clear;
global w
i=1;
for w=100:10:1000;
options=odeset('relTol',1e-5,'Maxstep',1e-2);
    [T,Y]= ode45(@rigid,[0 10],[0 1 0 1 0 1 0 1],options)
    a(i)=Y(end,1);
    b(i)=Y(end,3);
    c(i)=Y(end,5);
    d(i)=Y(end,7);
    A1(i)=sqrt(a(i)^2+b(i)^2);
    A2(i)=sqrt(c(i)^2+d(i)^2);
    i=i+1;
    i
end
plot([100:10:1000],A1)