MATLAB实现弹道仿真

syms  y0 y1 y2 y3 y4 y5 x x0 x1 x2 x3 x4 x5 a Ma
t0=（x-x1）*（x-x2）*（x-x3）*（x-x4）*y0/（（x0-x1）*（x0-x2）*（x0-x3）*（x0-x4））   %对于阻力因数的线性差值，采用4阶拉格朗日插值法
t1=（x-x0）*（x-x2）*（x-x3）*（x-x4）*y1/（（x1-x0）*（x1-x2）*（x1-x3）*（x1-x4））
t2=（x-x1）*（x-x0）*（x-x3）*（x-x4）*y2/（（x2-x1）*（x2-x0）*（x2-x3）*（x2-x4））
t3=（x-x1）*（x-x2）*（x-x0）*（x-x4）*y3/（（x3-x1）*（x3-x2）*（x3-x0）*（x3-x4））
t4=（x-x1）*（x-x2）*（x-x3）*（x-x0）*y4/（（x4-x1）*（x4-x2）*（x4-x3）*（x4-x0））
Cx0=t0+t1+t2+t3+t4
Cx15=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{a2468100.0430.05110.06510.08470.112}）
Cx21=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{a2468100.0360.04360.05580.07360.0973}）
Cx27=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{a2468100.03080.03720.04810.06410.0849}）
Cx33=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{a2468100.02650.03230.04190.0560.0746}）
Cx40=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{a2468100.02220.02720.03560.04780.0644}）
Cx=subs（Cx0{xx0x1x2x3x4y0y1y2y3y4}{Ma1.52.12.73.34.0Cx15Cx21Cx27Cx33Cx40}）
f0=（x-x1）*（x-x2）*（x-x3）*（x-x4）*（x-x5）*y0/（（x0-x1）*（x0-x2）*（x0-x3）*（x0-x4）*（x0-x5））   %对于升力因数的线性差值，采用5阶拉格朗日插值法
f1=（x-x0）*（x-x2）*（x-x3）*（x-x4）*（x-x5）*y1/（（x1-x0）*（x1-x2）*（x1-x3）*（x1-x4）*（x1-x5））
f2=（x-x1）*（x-x0）*（x-x3）*（x-x4）*（x-x5）*y2/（（x2-x1）*（x2-x0）*（x2-x3）*（x2-x4）*（x2-x5））
f3=（x-x1）*（x-x2）*（x-x0）*（x-x4）*（x-x5）*y3/（（x3-x1）*（x3-x2）*（x3-x0）*（x3-x4）*（x3-x5））
f4=（x-x1）*（x-x2）*（x-x3）*（x-x0）*（x-x5）*y4/（（x4-x1）*（x4-x2）*（x4-x3）*（x4-x0）*（x4-x5））
f5=（x-x1）*（x-x2）*（x-x3）*（x-x4）*（x-x0）*y5/（（x5-x1）*（x5-x2）*（x5-x3）*（x5-x4）*（x5-x0））
Cya0=f0+f1+f2+f3+f4+f5
Cya15=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.03020.03040.03060.03090.03110.0313}）
Cya20=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.02790.02800.02840.02860.02880.029}）
Cya25=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.02610.02640.02670.02690.02720.0274}）
Cya30=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.02470.02480.02510.02540.02570.0259}）
Cya35=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.02260.02270.02310.02330.02360.0238}）
Cya40=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{a12468100.02090.0210.02130.02160.02190.0221}）
Cya=subs（Cya0{xx0x1x2x3x4x5y0y1y2y3y4y5}{Ma1.522.533.54Cya15Cya20Cya25Cya30Cya35Cya40}）
mu=zeros（1100）;           %为迭代量预置空间，
Y=zeros（1100）;
theta