simulate an object thrown in 3d with drag force using ode45 in matlab

main.m
dx/dt = vx
dy/dt = vy
dz/dt = vz
dvx/dt = -(F_d / m) * vx
dvy/dt = -g - (F_d / m) * vy
dvz/dt = -(F_d / m) * vz

where:
F_d = 0.5 * rho * v^2 * C_d * A
rho = air density
v = velocity
C_d = drag coefficient
A = cross-sectional area
m = object mass
g = acceleration due to gravity
279 chars
16 lines

main.m
function dYdt = obj_motion(t, Y, rho, C_d, A, m)
    % unpack Y
    x = Y(1);
    y = Y(2);
    z = Y(3);
    vx = Y(4);
    vy = Y(5);
    vz = Y(6);
    
    % constants
    g = 9.81; % m/s^2
    
    % calculate drag force
    v = sqrt(vx^2 + vy^2 + vz^2);
    F_d = 0.5 * rho * v^2 * C_d * A;
    
    % equations of motion
    dxdt = vx;
    dydt = vy;
    dzdt = vz;
    dvxdt = -(F_d / m) * vx;
    dvydt = -g - (F_d / m) * vy;
    dvzdt = -(F_d / m) * vz;
    
    % pack up output
    dYdt = [dxdt; dydt; dzdt; dvxdt; dvydt; dvzdt];
end
546 chars
28 lines

main.m
% initial conditions
Y0 = [0; 0; 0; v0*sin(theta)*cos(phi); v0*sin(theta)*sin(phi); v0*cos(theta)];

% constants
rho = 1.225; % kg/m^3
C_d = 0.47;
A = 0.01; % m^2
m = 0.1; % kg

% time span
tspan = [0 10];

% call ode45 to solve the equations
[t, Y] = ode45(@(t, Y) obj_motion(t, Y, rho, C_d, A, m), tspan, Y0);
312 chars
15 lines