function xdot = STOL_EOM(t,x)
global e1 e2 e3 rho m g S b c AR WE Power ...
GeneralizedInertia GeneralizedInertia_Inv ...
CD0 e CL_Trim CLalpha CLq CLde Cmde ...
CYbeta CYp CYr Cmalpha Cmq Clbeta Clp Clr Cnbeta Cnp Cnr ...
CYda CYdr Clda Cldr Cnda Cndr ...
P_dynamic = (1/2)*rho*norm(V)^2;
RIB = expm(psi*hat(e3))*expm(theta*hat(e2))*expm(phi*hat(e1));
LIB = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
CL = CL_Trim + CLalpha*alpha + CLq*((q*c)/(2*norm(V))) + CLde*de;
CD = CD0 + (CL^2)/(e*pi*AR);
temp = expm(-alpha*hat(e2))*expm(beta*hat(e3))*[-CD; 0; -CL];
CY = CYbeta*beta + CYp*((b*p)/(2*norm(V))) + CYr*((b*r)/(2*norm(V))) + ...
X = P_dynamic*S*CX + Thrust;
Cl = Clbeta*beta + Clp*((b*p)/(2*norm(V))) + Clr*((b*r)/(2*norm(V))) + ...
Cm = Cmalpha*alpha + Cmq*((c*q)/(2*norm(V))) + Cmde*de;
Cn = Cnbeta*beta + Cnp*((b*p)/(2*norm(V))) + Cnr*((b*r)/(2*norm(V))) + ...
temp = GeneralizedInertia*[VE; omega];
LinearMomentum = temp(1:3);
AngularMomentum = temp(4:6);
RHS = [cross(LinearMomentum,omega) + Force; ...
cross(AngularMomentum,omega) + Moment];
temp = GeneralizedInertia_Inv*RHS;
VDot = VEDot + cross(omega,RIB'*WE);
xdot = [XDot; ThetaDot; VDot; omegaDot];
====================================script=================
global e1 e2 e3 rho m g S b c AR Inertia
P_dynamic = (1/2)*rho*u0^2;
CL_Trim = W/(P_dynamic*S);
CD0 = CD_Trim - (CL_Trim^2)/(e*pi*AR);
Thrust_Trim = P_dynamic*S*CD_Trim;
y0 = [X0; Theta0; V0; omega0];
[t,y] = ode45('STOL_EOM',[0:0.1:t_final]',y0);
ylabel('Angular Velocity')