请大家看看我的程序 我搞不明白哪里不对啊

jev

  运行起来会出错 请高手看看
function  dy=srotor(t,x)
global r
R=100;B=10;L=20;CR=2;M=10;K=10^6;e=0;g=10^3;I=10^(-3);q=0.1;ps=5*10^5;pb=10^5;w=400;us=0.5;u0=200;pc=0;
q1=R/(R+B);q2=B/CR;q3=sqrt(ps)/(g*w*R);q4=R/L;q5=R/CR;q6=g*L*w*R/I;q7=L/CR;q8=R/B;q9=B/L;q10=K/(M*w);q11=(ps*L*R)/(M*CR*w^2);q12=e/CR;
dy=zeros(10,1);
dy(1)=x(3);  
dy(2)=x(4);     
dy(3)=-q10*x(1)-3.14*q11*pc+q12*cost;
dy(4)=-q10*x(2)-3.14*q11*ps+q12*sint;
dy(5)=(1/(x(1)^2+x(2)^2-2))*((q1*(q3)^2)*(x(2)*x(12)-x(1)*x(13))+q1*(x(2)*x(6)-x(2)*x(7))+q2*(x(7)*x(9)-x(6)*x(10))+2*(q3*q4*(q3*q4*(1+0.24*x(1)+0.971*x(2))*sqrt(2*(1-x(11)+0.24*x(12)+0.971*x(13))/(1+q))*(x(5)-0.24*x(6)-0.971*x(7)-us)+q3*q4*(1-0.24*x(1)-0.971*x(2))*sqrt(2*(1-x(11)-0.24*x(12)-0.971*x(13))/(1+q))*(x(5)+0.24*x(6)+0.971*x(7)-us))-q2*(0.142*q8*(x(8)-0.24*x(9)-0.971*x(10)-x(5)+0.24*x(6)+0.971*x(7))*sqrt((x(8)-0.24*x(9)-0.971*x(10)-x(5)+0.24*x(6)+0.971*x(7))^2+(1.26*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2))/(1+q)+0.142*q8*(x(8)+0.24*x(9)+0.971*x(10)-x(5)-0.24*x(6)-0.971*x(7))*sqrt((x(8)+0.24*x(9)+0.971*x(10)-x(5)-0.24*x(6)-0.971*x(7))^2+(1.26*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2))/(1+q))+(0.0395*q5(x(5)-0.24*x(6)-0.971*x(7))*sqrt(((2*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2)/(1+q))*(x(5)-0.24*x(6)-0.971*x(7)))*(((2*q6*(1+0.24*x(1)+0.971*x(2)))/(1+0.24*x(1)+0.971*x(2)+q7)))*sqrt(((2*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2)/(1+q))*(x(5)-0.24*x(6)-0.971*x(7))))^(-0.25)+0.0395*q5(x(5)+0.24*x(6)+0.971*x(7))*sqrt(((2*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2)/(1+q))*(x(5)+0.24*x(6)+0.971*x(7)))*(((2*q6*(1-0.24*x(1)-0.971*x(2)))/(1-0.24*x(1)-0.971*x(2)+q7)))*sqrt(((2*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2)/(1+q))*(x(5)+0.24*x(6)+0.971*x(7)))^(-0.25))+x(1)*T1c+y*T1s);
dy(6)=(1/(x(1)^2+x(2)^2-2))*(x(1)*((q1*(q3)^2)*(x(2)*x(12)-x(1)*x(13))+q1*(x(2)*x(6)-x(1)*x(7))+q2*(x(7)*x(9)-x(6)*x(10))+2*(q3*q4*(D1+D2)-q2*(A1+A2)+(F1+F2)))+(2-x(2)^2)*(x(13)*q1*(q3)^2+q1*x(5)*x(7)+q2*(u0*x(10)-x(5)*x(10))+0.48*(q3*q4*(D2-D1)-q2*(A2-A1)+(F2-F1)))+x(1)*x(2)*(-1*x(12)*q1*(q2)^2+q1*x(5)*x(6)+q2*(x(5)*x(9)-u0*x(9))+1.942*(q3*q4*(D2-D1)-q2*(A2-A1)+(F2-F1))));
dy(7)=(1/(x(1)^2+x(2)^2-2))*(x(1)*T1+(2-x(1)^2)*T1s+x(1)*x(2)*T1c);
dy(8)=-1/T2;
dy(9)=-(2*(A1+A2-(2*q9+1)*(B1+B2)));
dy(10)=-(-1*x(12)*(q3)^2+u0*x(6)-2*x(8)*x(9)+1.942*(A2-A1+(2*q9+1)*(B2-B1)));

A1=0.142*q8*(x(8)-0.24*x(9)-0.971*x(10)-x(5)+0.24*x(6)+0.971*x(7))*sqrt((x(8)-0.24*x(9)-0.971*x(10)-x(5)+0.24*x(6)+0.971*x(7))^2+(1.26*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2))/(1+q);
A2=0.142*q8*(x(8)+0.24*x(9)+0.971*x(10)-x(5)-0.24*x(6)-0.971*x(7))*sqrt((x(8)+0.24*x(9)+0.971*x(10)-x(5)-0.24*x(6)-0.971*x(7))^2+(1.26*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2))/(1+q);
B1=0.0395*q8*(1-x(8)+0.24*x(9)+0.971*x(10))*sqrt((1-x(8)+0.24*x(9)+0.971*x(10))^2+(0.0845*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2)/(1+q))*(((2*q6*q9)/(1+q9))*sqrt((1-x(8)+0.24*x(9)+0.971*x(10))^2+(0.0845*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2)/(1+q)))^(-0.25);
B2=0.0395*q8*(1-x(8)-0.24*x(9)-0.971*x(10))*sqrt((1-x(8)-0.24*x(9)-0.971*x(10))^2+(0.0845*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2)/(1+q))*(((2*q6*q9)/(1+q9))*sqrt((1-x(8)-0.24*x(9)-0.971*x(10))^2+(0.0845*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2)/(1+q)))^(-0.25);
D1=q3*q4*(1+0.24*x(1)+0.971*x(2))*sqrt(2*(1-x(11)+0.24*x(12)+0.971*x(13))/(1+q))*(x(5)-0.24*x(6)-0.971*x(7)-us);
D2=q3*q4*(1-0.24*x(1)-0.971*x(2))*sqrt(2*(1-x(11)-0.24*x(12)-0.971*x(13))/(1+q))*(x(5)+0.24*x(6)+0.971*x(7)-us);
F1=0.0395*q5*(x(5)-0.24*x(6)-0.971*x(7))*sqrt((2*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2/(1+q))*(x(5)-0.24*x(6)-0.971*x(7))^2)*((2*q6*(1+0.24*x(1)+0.971*x(2))/(1+0.24*x(1)+0.971*x(2)+q7))*sqrt((2*(1-x(11)+0.24*x(12)+0.971*x(13))*(q3)^2/(1+q))*(x(5)-0.24*x(6)-0.971*x(7))^2))^(-0.25);
F2=0.0395*q5*(x(5)+0.24*x(6)+0.971*x(7))*sqrt((2*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2/(1+q))*(x(5)+0.24*x(6)+0.971*x(7))^2)*((2*q6*(1-0.24*x(1)-0.971*x(2))/(1-0.24*x(1)-0.971*x(2)+q7))*sqrt((2*(1-x(11)-0.24*x(12)-0.971*x(13))*(q3)^2/(1+q))*(x(5)+0.24*x(6)+0.971*x(7))^2))^(-0.25);
G1=q3*q4*(1+0.24*x(1)+0.971*x(2))*(sqrt(2*(1-x(11)+0.24*x(12)+0.971*x(13))/(1+q))-sqrt(2*(x(11)-0.24*x(12)-0.971*x(13)-pb)/(1+q)));
G2=q3*q4*(1-0.24*x(1)-0.971*x(2))*(sqrt(2*(1-x(11)-0.24*x(12)-0.971*x(13))/(1+q))-sqrt(2*(x(11)+0.24*x(12)+0.971*x(13)-pb)/(1+q)));
T1=(q1*(q3)^2)*(x(2)*x(12)-x(1)*x(13))+q1*(x(2)*x(6)-x(1)*x(7))+q2*(x(7)*x(9)-x(6)*x(10))+2*(q3*q4*(D1+D2)-q2*(A1+A2)+(F1+F2));
T1c=x(13)*q1*(q3)^2+q1*x(5)*x(7)+q2*(u0*x(10)-x(5)*x(10))+0.48*(q3*q4*(D2-D1)-q2*(A2-A1)+(F2-F1));
T1s=-1*x(12)*q1*(q2)^2+q1*x(5)*x(6)+q2*(x(5)*x(9)-u0*x(9))+1.942*(q3*q4*(D2-D1)-q2*(A2-A1)+(F2-F1));
T2=2*(A1+A2-(2*q9+1)*(B1+B2));
T2c=x(13)*(q3)^2+2*x(8)*x(10)-u0*x(10)+0.48*(A2-A1+(2*q9+1)*(B2-B1));
T2s=-1*x(12)*(q3)^2+u0*x(9)-2*x(8)*x(9)+1.942*(A2-A1+(2*q9+1)*(B2-B1));
T3=2*(G1+G2);
T3c=-1*q1*x(5)*x(2)+q1*x(7)+q2*x(10)+0.48*(G2-G1);
T3s=q1*x(5)*x-q1*x(6)-q2*x(9)+1.942*(G2-G1);
[t,x]=ode45('srotor',[0 5],[0 ;0 ;0 ;0 ;0 ;0 ;0 ;0 ;0 ;0 ;0.5 ;0.5 ;0])