Вывод уравнений движения ЛА типа "самолет"
Допущения:
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={vk[t],0,0}] (*Вектор земной скорости в траекторной системе координат*)

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

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

ukpd1=D[R[[1]],t]==vg[[1]]

ukpd2=D[R[[2]],t]==vg[[2]]

ukpd3=D[R[[3]],t]==vg[[3]]

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

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

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

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

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

udpd1=(massa*(D[v,t]+Cross[w,v]))[[1]]==f[[1]]

udpd2=(massa*(D[v,t]+Cross[w,v]))[[2]]==f[[2]]

udpd3=(massa*(D[v,t]+Cross[w,v]))[[3]]==f[[3]]

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

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

ukvd1=w[[1]]==proekc[[1]]

ukvd2=w[[2]]==proekc[[2]]

ukvd3=w[[3]]==proekc[[3]]

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

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] (*Полная производная по времени главного вектора момента количеств движения*)

Полная производная главного вектора момента количеств движения равна главному моменту внешних сил, действующих на тело.

udvd1=lchudvd[[1]]==M[[1]]

udvd2=lchudvd[[2]]==M[[2]]

udvd3=lchudvd[[3]]==M[[3]]

MatrixForm[polnsyst={ukpd1,ukpd2,ukpd3,udpd1,udpd2,udpd3,ukvd1,ukvd2,ukvd3,udvd1,udvd2,udvd3}]

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

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

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

VectorSostOpor={L0[t],H0[t],dZ0[t],vx0[t],vy0[t],vz0[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],vx[t]->vx0[t],vy[t]->vy0[t],vz[t]->vz0[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],Δvx[t],Δvy[t],Δvz[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],vx[t]->vx0[t]+Δvx[t],vy[t]->vy0[t]+Δvy[t],vz[t]->vz0[t]+Δvz[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,14} ]]
(*dfi[t]/dxi*)