The drag force should always be in the opposite direction of current velocity direction. You don't have that. You always have drag in one direction that is not dependent on the current velocity direction. E.g.,
athrust(i) = (fthrust(i) - fdrag(i - 1)) ./ mtotal(i);
:
fdrag(i) = ((vel(i - 1)).^2) * .0075;
Looking at the above code, there is nothing to change the sign of the force due to drag based on the current direction of velocity. You need to add code to change that. E.g., maybe something like this:
fdrag(i) = vel(i-1) * abs(vel(i-1)) * .0075;
EDIT:
If you want the drag force plot to match the "correct" one for sign, then incorporate the sign into the force and flip the sign on the application. E.g., flip the sign on fdrag in both places:
athrust(i) = (fthrust(i) + fdrag(i - 1)) ./ mtotal(i);
:
fdrag(i) = -vel(i-1) * abs(vel(i-1)) * .0075;
Best Answer