MATLAB: I’m having trouble with the drag on the rocket projectile.

dragMATLABplotrocket

The graphs aren't coming out right once I add drag. Any help would be appreciated. (Attached is my plot of drag and the plot I am trying for.)
Edit: Added abs(vel(i-1)
        % PROGRAM DESCRIPTION
        % Calculates the trajectory of a rocket's flight with no air resistance. The
        % program will then plot several graphs including acceleration vs time,
        % force of thrust vs time, and total energy vs time.
        % Initial Conditions
        % When time is 0, the elevation of the rocket is 1 meter and the mass of the rocket is 150kg. There is a time increment of .01 and acceleration
        % due to gravity is constant.
        % PROGRAM VARIABLES
        % ag       Acceleration due to gravity (m/s^2)
        % athrust  Acceleration due to thrust
        % atotal   Total acceleration (m/s^2)
        % dtime    Length of time step
        % fdrag    Force of drag (N)
        % fthrust  Force of thrust (N)
        % height   Elevation above ground (m)
        % ke       Kinetic Energy (J)
        % mfuel    Mass of fuel (kg)
        % mtotal   Total mass of rocket (kg)
        % pe       Potential energy (J)
        % time     Elapsed time (s)
        % te       Total energy (J)
        % vel      Velocity (m/s)
        % weight   Weight of rocket (kg)
        clear all
        close all
        time(1) = 0;
        dtime = .01;
        fthrust(1) = 1560;
        height(1) = 1;
        weight(1) = 150;
        mfuel(1) = 132;
        mrocket = 18;
        vel(1) = 0; 
        ag = -9.81;
        atotal(1) = 0;
        ftotal = 0;
        fdrag = 0;
        i = 1;
        % Calculates the variables involved in the launch of the rocket. Loop will
        % run until rocket reaches the ground.
        while height(i) > 0
            i = i + 1;
            time(i) = time(i-1) + dtime;
            % If the rocket still has fuel, the rocket will have a thrust of 1560N.
            % Otherwise, there is no fuel nor thrust.
            if (mfuel(i - 1) - 5 * dtime) > 0 
                mfuel(i) = mfuel(i - 1) - 5 * dtime;
                fthrust(i) = fthrust(1);
            else
                mfuel(i) = 0;
                fthrust(i) = 0;
            end
            mtotal(i) = mrocket + mfuel(i);
            athrust(i) = (fthrust(i) - fdrag(i - 1)) ./ mtotal(i);
            % If there is not any fuel in the rocket, only gravity will accelerate
            % it. With fuel, the acceleration includes the acceleration due to the
            % thrusters.
            if mfuel(i) == 0
                atotal(i) = ag;
            else
                atotal(i) = (athrust(i - 1) + ag);
            end
            fdrag(i) = vel(i - 1) * abs(vel(i - 1)) * .0075;
            vel(i) = vel(i - 1) + atotal(i) * dtime;
            height(i) = height(i - 1) + vel(i - 1) * dtime + (.5 * atotal(i - 1) * dtime .^ 
        2);
            ke(i) = .5 * (mtotal(i) * (vel(i - 1) .^ 2));
            pe(i) = -mtotal(i - 1) * ag * height(i - 1);
            te(i) = ke(i) + pe(i);
        end
        %Plots Height, Velocity, and Acceleration vs Time
        subplot(2,2,1);
        grid; 
        plot(time, atotal, 'k');
        title('Acceleration: With Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Acceleration (m/s^2)');
        subplot(2,2,2);
        grid; 
        plot(time, vel, 'k');
        title('Velocity: With Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Velocity (m/s)');
        subplot(2,2,3);
        grid; 
        plot(time, height, 'k');
        title('Altitude: With Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Altitude (m)');
        subplot(2,2,4);
        grid; 
        plot(time, mtotal, 'k');
        title('Total Mass: With Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Mass (kg)');
        %Plots Thrust Force vs Time in new figure
        figure; 
        plot(time, fthrust, 'r');
        grid; 
        title('Force of Thrust: With Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Thrust (N)'); 
        %Plots Energy vs Time in new figure
        figure;
        plot(time, ke);
        hold on;
        plot(time, pe, 'g');
        plot(time, te, 'r');
        grid; 
        hold off;
        title('Energy: Kinetic, Potential, Total');
        xlabel('Time (seconds)');
        ylabel('Energy (J)');
        legend('Kinetic', 'Potential', 'Total');
        % Plots Drag Force vs Time in new figure
        figure;
        plot(time, fdrag);
        grid;
        title('Force of Drag due to Wind Resistance');
        xlabel('Time (seconds)');
        ylabel('Drag Force (N)');

Best Answer

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;
Related Question