function dudt = f_withdrag(t, u)
% equations of motion for baseball in flight with turbulent air resistance
  
    rho_air  = 1.196;   % kg/m^3, density of dry air, 21 C, sea level   (Fenway)
    C = 0.3;            % drag coefficient for baseball (rough sphere)
    g = 9.81;           % acceleration due to gravity in m/s^2
    r = 0.0375;         % radius of baseball in m (3.75 cm)
    A = pi*r^2;         % cross-sectional area of baseball in m^2
    m = 0.145;          % mass of baseball in kg (145 gm
    mu = rho_air*C*A/2; % coefficient of nonlinear |v|^2 term, in mks units

    % extract components of vector u into variables x, y, vx, vy
    x  = u(1);
    y  = u(2);
    vx = u(3);
    vy = u(4);

    % compute rates of change of those components
    dxdt = vx;
    dydt = vy;
    dvxdt = -mu/m * vx * sqrt(vx^2 + vy^2);
    dvydt = -mu/m * vy * sqrt(vx^2 + vy^2) - g;

    % pack rates of change back into vector dudt, return value of function    
    dudt = [dxdt; dydt; dvxdt; dvydt];
end