Code:

```
S = 0.017; % Reference Area, m^2
AR = 0.86; % Wing Aspect Ratio
e = 0.9; % Oswald Efficiency Factor;
m = 0.003; % Mass, kg
g = 9.8; % Gravitational acceleration, m/s^2
rho = 1.225; % Air density at Sea Level, kg/m^3
CLa = 3.141592 * AR/(1 + sqrt(1 + (AR / 2)^2));
% Lift-Coefficient Slope, per rad
CDo = 0.02; % Zero-Lift Drag Coefficient
epsilon = 1 / (3.141592 * e * AR);% Induced Drag Factor
CL = sqrt(CDo / epsilon); % CL for Maximum Lift/Drag Ratio
CD = CDo + epsilon * CL^2; % Corresponding CD
LDmax = CL / CD; % Maximum Lift/Drag Ratio
Gam = -atan(1 / LDmax); % Corresponding Flight Path Angle, rad
V = sqrt(2 * m * g /(rho * S * (CL * cos(Gam) - CD * sin(Gam))));
% Corresponding Velocity, m/s
Alpha = CL / CLa; % Corresponding Angle of Attack, rad
% Oscillating Glide due to Zero Initial Flight Path Angle
xo = [V;0;H;R];
[tb,xb] = ode23('EqMotion',tspan,xo);
function xdot = EqMotion(t,x)
% Fourth-Order Equations of Aircraft Motion
global CL CD S m g rho
V = x(1);
Gam = x(2);
q = 0.5 * rho * V^2; % Dynamic Pressure, N/m^2
xdot = [(-CD * q * S - m * g * sin(Gam)) / m
(CL * q * S - m * g * cos(Gam)) / (m * V)
V * sin(Gam)
V * cos(Gam)];
```

I've also learned how to solve differential equations using the Runge-Kutta Fourth Order Method. However, when I try to solve for the next timestep using RK4, the values do not correspond with the graphs in MATLAB, which means I'm doing something wrong.

So my questions are, first, am I correct in assuming that the four lines of xdot correspond to velocity, gamma, x position, and y position respectively? If this is true, can I use the velocity and gamma equations to solve for v(n+1) and g(n+1) using the Runge-Kutta method?

I'm sorry for my lack of knowledge, but I hope you can help in answering my question. Please feel free to ask any questions and I'll do my best to clarify.