Simple Drag Simulation

nstep = 10000;
dt = 0.01; % centisecond step
% m = [1.15 1.0]; % mass, kg
% m = [2883 2394]; % Yak-1b loaded/unloaded masses (Wikipedia)
m = [2656.4 2362.6]; % Estimated weights

v = zeros(nstep,length(m)); % Velocity
D = zeros(nstep,length(m)); % Drag force

% rho = 1.2754; % Density of air @ STP (IUPAC), kg/m^3
% rho = 1.315; % Air density @ -15C, 300m, 50% RH
rho = 1.07; % Air density @ -15C, 2000m, 50% RH
% cd = 0.027; % Drag coefficient of Cessna 172
cd = 0.025; % Hypothetical WW2 fighter
% A = 17.2; % Wing area, m^2
A = 17.2; % Yak-1b (Wikipedia)
drag = @(x) 0.5 * rho * cd * A * x.^2;

vmax = 550 * 1000 / 3600;

v(1,:) = 650 * 1000 / 3600;
D(1,:) = drag(v(1));
breakstep = nstep;
for i=2:nstep
    for j=1:length(m)
        v(i,j) = max(v(i-1,j) - D(i-1,j) / m(j) * dt,vmax);
        D(i,j) = drag(v(i,j));
    end
    if all(v(i,:)<=vmax)
        breakstep = i;
        break;
    end
end
v(breakstep+1:end,:) = [];
D(breakstep+1:end,:) = [];