Вывод уравнений движения ЛА типа "самолет"
Допущения:
1) Конструкция ЛА абсолютно жесткая (ЛА - твердое тело);
2) Атмосфера спокойная (W=0 - скорость ветра);
3) Масса ЛА постоянная т.к. заданное время переходных процессов меньше 60 с;
4) Земля плоская и не вращается;
5) Угол тангажа не равен +/-90 градусов;
6) Спарведлива гипотеза стационарности (согласно этой гипотезе, аэродинамическая сила и аэродинамический момент в нестационарном режиме полета полностью определяются величинами параметров в данный момент времени).
Системы координат и матрицы перехода

MatrixForm[M1={{Cos[Θ[t]],-Sin[Θ[t]],0},{Sin[Θ[t]],Cos[Θ[t]],0},{0,0,1}}] (*Матрица поворота на угол наклона траектории*)

MatrixForm[M2={{Cos[Ψ[t]],0,Sin[Ψ[t]]},{0,1,0},{-Sin[Ψ[t]],0,Cos[Ψ[t]]}}] (*Матрица поворота на угол пути*)

MatrixForm[TN=M2.M1] (*Матрица перехода от траекторной системы координат к нормальной системе координат *)

MatrixForm[ReplaceAll[ TN,{Φ[t]->0,Θ[t]->0}]]  (*Проверка подстановкой нулевых углов поворотов. Должна получиться единичная матрица*)

MatrixForm[MM1={{Cos[ψ[t]],0,Sin[ψ[t]]},{0,1,0},{-Sin[ψ[t]],0,Cos[ψ[t]]}}] (*Матрица поворота на угол рыскания*)
MatrixForm[Transpose[MM1]]

MatrixForm[MM2={{Cos[ϑ[t]],-Sin[ϑ[t]],0},{Sin[ϑ[t]],Cos[ϑ[t]],0},{0,0,1}}](*Матрица поворота на угол тангажа*)
MatrixForm[Transpose[MM2]]

MatrixForm[MM3={{1,0,0},{0,Cos[γ[t]],-Sin[γ[t]]},{0,Sin[γ[t]],Cos[γ[t]]}}] (*Матрица поворота на угол крена*)
MatrixForm[Transpose[MM3]]

MatrixForm[SN=MM1.MM2.MM3] (*Матрица перехода от связанной системы координат к нормальной системе координат*)

MatrixForm[ReplaceAll[ SN,{ϑ[t]->0,ψ[t]->0,γ[t]->0} ]]  (*Проверка подстановкой нулевых углов поворотов. Должна получиться единичная матрица*)

MatrixForm[SM1={{Cos[α[t]],-Sin[α[t]],0},{Sin[α[t]],Cos[α[t]],0},{0,0,1}}] (*Матрица поворота на угол атаки*)

MatrixForm[SM2={{Cos[β[t]],0,Sin[β[t]]},{0,1,0},{-Sin[β[t]],0,Cos[β[t]]}}] (*Матрица поворота на угол скольжения*)

MatrixForm[SS=Transpose[SM2.SM1]] (*Матрица перехода от скоростной системы координат к связанной системе координат*)

MatrixForm[ReplaceAll[ SS,{α[t]->0,β[t]->0}]]  (*Проверка подстановкой нулевых углов поворотов. Должна получиться единичная матрица*)

Кинематические уравнения поступательного движения

R={L[t],H[t],dZ[t]}(*Радиус вектор центра масс ЛА, L[t] - дальность полета, H[t] - высота полета, dZ[t] - боковое отклонение*)

vg={vgx[t],vdy[t],vgz[t]} (*Вектор земной сокрости центра масс ЛА в нормальной системе координат*)

MatrixForm[D[R,t]]==MatrixForm[vg] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

MatrixForm[vzem={V[t],0,0}] (*Вектор земной скорости в траекторной системе координат*)

MatrixForm[vg=TN.vzem] (*Проецируем вектор земной скорости на нормальные оси координат*)

MatrixForm[udpdk1=Table[ D[R,t] [[i]]==vg[[i]],{i,1,3}]]
 (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

(*MatrixForm[udpdk2=Table[vzem[[i]]==(SS.vzem)[[i]],{i,1,3}]]*)

Динамические уравнения поступательного движения

v={vx[t],vy[t],vz[t]} (*Вектор воздушной скорости ЛА в связанной системе координат*)

f={fx,fy,fz} (*Вектор результирующей сил, действующих на ЛА в полете, в связанной системе координат*)

w={wx[t],wy[t],wz[t]} (*Вектор угловой скорости ЛА в связанной системе координат*)

MatrixForm[udpdd=Table[(massa*(D[v,t]+Cross[w,v]))[[i]]==f[[i]],{i,1,3}]]
 (*Динамическое уравнение движения центра масс самолета в векторной форме*)

MatrixForm[udpdk2=Table[v[[i]]==(SS.vzem)[[i]],{i,1,3}]]

Кинематические уравнения вращательного движения

MatrixForm[proekc={D[γ[t],t],0,0}+Transpose[MM2.MM3].{0,D[ψ[t],t],0}+  Transpose[MM3].{0,0,D[ϑ[t],t]}]
(*Проецируем угловые скорости рыскания, тангажа и крена на оси связанной системы 
координат*)

MatrixForm[udvdk1=Table[w[[i]]==proekc[[i]],{i,1,3}]]

MatrixForm[udvdk2={Θ[t]==ϑ[t]-α[t],Ψ[t]==ψ[t]-β[t]}]

Динамические уравнения вращательного движения

MatrixForm[TI={{Ixx,-Ixy,-Ixz},{-Iyx,Iyy,-Iyz},{-Izx,-Izy,Izz}}] (*Тензор инерции*)

Ixz=Izx=Iyz=Izy=0; (*Ось OZ перпендикулярна плоскости симметрии XOY*)

Iyx=Ixy;

MatrixForm[TI]

MatrixForm[K=TI.w] (*Главный момент количеств движения в связанной системе координат*)

M={Mx,My,Mz} (*Результирующий момент аэродинамических сил и тяги в связанной системе координат(Момент силы тяжести, центробежных и кориолисовых сил инерции вокруг центра масс ЛА равен нулю)*)

lchudvd=D[K,t]+Cross[w,K] ;(*Полная производная по времени главного вектора момента количеств движения*)

MatrixForm[udvdd=Table[lchudvd[[i]]==M[[i]],{i,1,3}]]

ПОЛНЫЕ НЕЛИНЕЙНЫЕ УРАВНЕНИЯ ДВИЖЕНИЯ

MatrixForm[polnsyst1=Flatten[{udpdk1,udpdk2,udpdd,udvdk1,udvdk2,udvdd}]]

MatrixForm[polnsyst2=ReplaceAll[polnsyst1,
{Θ[t]->-α[t]+ϑ[t],Ψ[t]->-β[t]+ψ[t],vx[t]->udpdk2[[1,2]],vy[t]->udpdk2[[2,2]],vz[t]->udpdk2[[3,2]],
vx'[t]->D[udpdk2[[1,2]],t],vy'[t]->D[udpdk2[[2,2]],t],vz'[t]->D[udpdk2[[3,2]],t]}
 ]]
(*Заменяем проекции скорости на связанную СК проекциями на полусвязанную и заменяем углы наклона траектории и сноса*)

MatrixForm[polnsyst3=
{polnsyst2[[1]],polnsyst2[[2]],polnsyst2[[3]],polnsyst2[[7]],polnsyst2[[8]],polnsyst2[[9]],polnsyst2[[10]],polnsyst2[[11]],polnsyst2[[12]],polnsyst2[[15]],polnsyst2[[16]],polnsyst2[[17]]}]
(*Вырезаем лишние уравнения*)
polnsyst=polnsyst3;

Распишем вектор состояния системы

VectorSost=Flatten[{R,V[t],w,ψ[t],ϑ[t],γ[t],α[t],β[t]}] (*Вектор состояния системы*)

{L[t],H[t],dZ[t],V[t],wx[t],wy[t],wz[t],ψ[t],ϑ[t],γ[t],α[t],β[t]}

dVectorSost=D[VectorSost,t] (*Производная вектора состояния системы*)

VectorSostOpor={L0[t],H0[t],dZ0[t],V0[t],wx0[t],wy0[t],wz0[t],ψ0[t],ϑ0[t],γ0[t],α0[t],β0[t]} 
 (*Вектор состояния опорного движения системы*)
RMOpor={L[t]->L0[t],H[t]->H0[t],dZ[t]->dZ0[t],V[t]->V0[t],wx[t]->wx0[t],wy[t]->wy0[t],wz[t]->wz0[t],ψ[t]->ψ0[t],ϑ[t]->ϑ0[t],γ[t]->γ0[t],α[t]->α0[t],β[t]->β0[t]} ;

dVectorSostOpor=D[VectorSostOpor,t] (*Производная вектора состояния опорного движения системы*)

VectorSostVozm={ΔL[t],ΔH[t],ΔdZ[t],ΔV[t],Δwx[t],Δwy[t],Δwz[t],Δψ[t],Δϑ[t],Δγ[t],Δα[t],Δβ[t]}
 (*Вектор состояния возмущенного движения системы*)

dVectorSostVozm=D[VectorSostVozm,t] (*Производная вектора состояния возмущенного движения системы*)

RMPolnoe={L[t]->L0[t]+ΔL[t],H[t]->H0[t]+ΔH[t],dZ[t]->dZ0[t]+ΔdZ[t],V[t]->V0[t]+ΔV[t],wx[t]->wx0[t]+Δwx[t],wy[t]->wy0[t]+Δwy[t],wz[t]->wz0[t]+Δwz[t],ψ[t]->ψ0[t]+Δψ[t],ϑ[t]->ϑ0[t]+Δϑ[t],γ[t]->γ0[t]+Δγ[t],α[t]->α0[t]+Δα[t],β[t]->β0[t]+Δβ[t]}

MatrixForm[Koshi1=Flatten[Simplify[Solve[polnsyst,D[VectorSost,t],Sort->False]]]] (*Разрешаем систему относительно производных*)

MatrixForm[s1=Table[Koshi1[[i,2]],{i,1,12}]]
(*f(x,u)*)

MatrixForm[s2=ReplaceAll[s1,RMPolnoe]]
(*f(x0+Δx,u0+Δu)*)

MatrixForm[s3=ReplaceAll[s1,RMOpor]]
(*f(x0,u0)*)

MatrixForm[s4=Table[s2[[i]]-s3[[i]],{i,1,12} ]]
(*f(x0+Δx,u0+Δu)-f(x0,u0)*)

(*MatrixForm[s4=Table[FullSimplify[s4[[i]]], {i,1,12}]]*)
(*Упрощаем *)

MatrixForm[a=Table[D[ s4[[i]],VectorSostOpor[[j]] ],{i,1,12}, {j,1,12} ]]
(*dfi[t]/dxi*)