Projectile Motion with Air Resistance

In addition to a vertical gravitational force on a 2D projectile motion, there is also a friction force to a certain extent due to air resistance.
This frictional force is usually in the opposite direction to velocity and is proportional to the square of the velocity: $\vec{F}_r=-kv\vec{v}$ (Drag force, $F_D=-(1/2)c\rho A v^2\vec{v}/\vert\vec{v}\vert$ here, c is the drag coefficient, $\rho$ the air density, and A the projectile's cross-sectional area).

If we write Newton's 2 $^{nd}$ law as a vector in two dimensions,

$\displaystyle m\vec{a}$ $\displaystyle =$ $\displaystyle \vec{F}_{net}$

$\displaystyle m\frac{d^2\vec{r}}{dt^2}$ $\displaystyle =$ $\displaystyle m\vec{g}-kv\vec{v}$
and component wise (where $k/m=\gamma$ ):

$\displaystyle \frac{d^2x}{dt^2}= -\gamma\left( \sqrt{v_x^2+v_y^2}\right) v_x$ $\displaystyle ~~~\&~~~$ $\displaystyle \frac{dx}{dt}= v_x$

$\displaystyle \frac{d^2y}{dt^2}=-g-\gamma\left( \sqrt{v_x^2+v_y^2}\right) v_y$ $\displaystyle ~~~\&~~~$ $\displaystyle \frac{dy}{dt}= v_y$
Now, we have a set of equations.

$\displaystyle \frac{d^2x}{dt^2}= -\gamma\left( \sqrt{v_x^2+v_y^2}\right) v_x$	$\displaystyle ~~~\&~~~$	$\displaystyle \frac{dx}{dt}= v_x$
$\displaystyle \frac{d^2y}{dt^2}=-g-\gamma\left( \sqrt{v_x^2+v_y^2}\right) v_y$	$\displaystyle ~~~\&~~~$	$\displaystyle \frac{dy}{dt}= v_y$