用c语言龙格库塔解微分方程组,大神们谁会用龙格库塔法解微分方程组

我要用龙格库塔法解一个微分方程组，都是一阶的，不过有十五个方程，相互关联。用龙格库塔法的基本步骤是怎样的？

以下是我的步骤，大神们看看有什么错误。

先编写一个名为funct的m文件，内容如下：

function dy=funct(t,y)

dy=zeros(15,1);

dy(1)=(T-delG*sin(y(9))+0.5*den*y(13)^2*S*CXB0)/(m+lam11);

dy(2)=((G*(yG*sin(y(9))-xG*cos(y(9)))+0.5*den*y(13)^2*S*L*(CNa*y(14)+CNr*y(6)*L/y(13)+CNe*dele)-m*xG*y(1)*y(6))*(m*xG+lam26)-(0.5*den*y(13)^2*S*(CYa*y(14)+CYr*y(6)*L/y(13))-delG*cos(y(9))-m*y(1)*y(6))*(m*xG+lam66))/((m*xG+lam26)^2-(m+lam22)*(Jz+lam66));

dy(3)=((delG*cos(y(9))*sin(y(7))+0.5*den*y(13)^2*S*(CZb*y(15)+CZp*y(4)*L/y(13)+CZq*y(5)*L/y(13)+CZe*delr)+m*y(1)*y(5))*(Jy+lam55)+(0.5*den*y(13)^2*S*L*(CMb*y(15)+CMp*y(4)*L/y(13)+CMq*y(5)*L/y(13)+CMe*delr)-G*(zG*sin(y(9))+xG*cos(y(9))*sin(y(7)))+m*xG*y(1)*y(5))*(m*xG-lam35))/((m+lam33)*(Jy+lam55)-(m*xG-lam35)^2);

dy(4)=(CRp*0.5*den*y(13)^2*S*L+G*cos(y(9))*(yG*sin(y(7))+zG*cos(y(7)))+0.5*den*y(13)^2*S*L*(CRb*y(15)+CRp*y(4)*L/y(13)+CRq*y(5)*L/y(13)+CRe*delr+CRe*deld)+m*y(1)*(yG*y(5)+zG*y(6)))/(Jx+lam44);

dy(5)=((delG*cos(y(9))*sin(y(7))+0.5*den*y(13)^2*S*(CZb*y(15)+CZp*y(4)*L/y(13)+CZq*y(5)*L/y(13)+CZe*delr)+m*y(1)*y(5))*(m*xG-lam35)+(0.5*den*y(13)^2*S*L*(CMb*y(15)+CMp*y(4)*L/y(13)+CMq*y(5)*L/y(13)+CMe*delr)-G*(zG*sin(y(9))+xG*cos(y(9))*sin(y(7)))+m*xG*y(1)*y(5))*(m+lam33))/((m+lam33)*(Jy+lam55)-(m*xG-lam35)^2);

dy(6)=((G*(yG*sin(y(9))-xG*cos(y(9)))+0.5*den*VT^2*S*L*(CNa*y(14)+CNr*y(6)*L/VT+CNe*dele)-m*xG*y(1)*y(6))*(m+lam22)-(0.5*den*VT^2*S*(CYa*y(14)+CYr*y(6)*L/VT)-delG*cos(y(9))-m*y(1)*y(6))*(m*xG+lam66))/((m+lam22)*(Jz+lam66)-(m*xG+lam26)^2);

dy(7)=y(4)-(y(5)*cos(y(7))-y(6)*sin(y(7)))*tan(y(9));

dy(8)=(y(5)*cos(y(7))-y(6)*sin(y(7)))/cos(y(9));

dy(9)=y(5)*sin(y(7))+y(6)*cos(y(7));

dy(10)=y(1)*cos(y(9))*cos(y(8))+y(2)*(sin(y(8))*sin(y(7))-sin(y(9))*cos(y(8))*cos(y(7)))+y(3)*(sin(y(8))*cos(y(7))-sin(y(9))*cos(y(8))*sin(y(7)));

dy(11)=y(1)*sin(y(9))+y(2)*cos(y(9))*cos(y(7))-y(3)*cos(y(9))*sin(y(7));

dy(12)=y(2)*(cos(y(8))*sin(y(7))-sin(y(9))*sin(y(8))*cos(y(7)))+y(3)*(cos(y(8))*cos(y(7))-sin(y(9))*sin(y(8))*sin(y(7)))-y(1)*cos(y(9))*sin(y(8));

dy(13)=sqrt(y(1)^2+y(2)^2+y(3)^2);

dy(14)=atan(-y(2)/y(1));

dy(15)=asin(y(3)/y(13));

end

其中的CMp之类参数的我会定义好具体的数值。

接下来在命令窗口输入

[t,y]=ode45('funct',[0 0.1 20],[VT 0 0 0 0 0 0 0 0 0 0 0 0 0 0]);

plot(t,y(:,10));

结果却总是出错，有人能帮我吗？具体错在哪里？我的步骤有问题吗？