$$ m \ddot{x} = - b \dot{x}$$ leads to:
\begin{equation} v_x(t) = v_0 e^{-\frac{t}{\tau}} \\ x(t) = v_0 \tau \left(1 - e^{-t/\tau} \right) \end{equation}where $\tau = m/b$
$$m \ddot{y} = m g - b \dot{y}$$
This can be solved to obtain:
\begin{equation} v_y(t) = v_{\rm{ter}} + (v_{y0} - v_{\rm{ter}})e^{\frac{-t}{\tau}}\\ y(t) = v_{\rm{ter}} t + (v_{y0} - v_{\rm{ter}}) \tau \left( 1 - e^{\frac{-t}{\tau}}\right) \end{equation}where $v_{\rm{ter}} = mg/b$
Free Body Diagram for Linear Drag
$$m \ddot{y} = m g - b \dot{y}$$
When the two forces are equal, there is no more acceleration and we have reached terminal velocity: $v_{\rm{ter}} = mg / b$.
Our differential equation is then: \begin{equation} m \dot{v}_y = - b(v_y - v_{\rm{ter}}) \end{equation} Let $u = v_y - v_{\rm{ter}}$ then we can write: $$m \dot{u} = -b u$$ This familiar equation is solved with an exponential: $$u = A e^{-t/\tau}$$ with $\tau = m/b$
Replacing the $u$:
$$v_y - v_{\rm{ter}} = A e^{-t/\tau}$$At $t = 0$, the velocity is $v_{y0}$ so the constant can be found to be: $$A = v_y - v_{\rm{ter}}$$
Since $$ y(t) = \int_0^t v_y(t) dt $$ obtaining the vertical position is done by performing the integral of \eqref{eq:vytlinear}.
\begin{equation} y(t) = v_{\rm{ter}} t + (v_{y0}-v_{\rm{ter}})\tau \left(1-e^{-t/\tau} \right) \end{equation}
Two dimensions - Two separated equations.
\begin{align} x(t) & = v_0 \tau \left(1 - e^{-t/\tau} \right) \\ y(t) & = (v_{y0} + v_{\rm{ter}}) \tau \left( 1 - e^{\frac{-t}{\tau}}\right) - v_{\rm{ter}} t \end{align}(note: for the following the direction of $y$ is now positive-up. )
If we want to express $y(x)$, that is, eliminate time $t$ from the equations, we can obtain:
\begin{equation} y = \frac{v_{y0}+ v_{\rm{ter} } } {v_{x0} } x + v_{\rm{ter}} \tau \ln \left(1-\frac{x}{v_{x0} \tau} \right) \label{eq:positionvstime-linear-drag} \end{equation}How can we make sense of this? (Don't even bother trying to solve for $R$)
Since $\tau = m / b$ and $v_\rm{ter} = mg /b$, we can expect both of these to be large. Thus $\frac{R}{v_{x0} \tau}$ from this term: $$\ln \left(1-\frac{R}{v_{x0} \tau} \right)$$ is small.
Using \begin{equation} \ln (1 - \epsilon) = - \left( \epsilon + \frac{1}{2} \epsilon^2 + \frac{1}{3}\epsilon^3 + \ldots \right) \end{equation}
We can arrive at:
\begin{equation} \left[\frac{v_{y0}+ v_{\rm{ter} } } {v_{x0} }\right] R - v_{\rm{ter}} \tau \left[ \frac{R}{v_{x0} \tau} + \frac{1}{2} \left(\frac{R}{v_{x0} \tau} \right)^2 + \frac{1}{3} \left( \frac{R}{v_{x0} \tau} \right)^3 \right] = 0 \end{equation}This is a lot nicer to deal with. It can be rewritten as a 2nd order polynomial:
\begin{equation} R = \frac{2 v_{x0} v_{y0}}{g} - \frac{2}{3 v_{x0} \tau}R^2 \end{equation}(we've replaced $v_{\rm{ter} }\tau$ with $g$)
If $\tau$ is large, i.e. little to no damping, then the $R^2$ term can be ignored, and we have the familiar range equation from physics in a vacuum.
\begin{equation} R \approx \frac{2v_{x0} v_{y0}}{g} = R_{\rm{vac}} \end{equation}In cartesian coordinates: $$ \begin{align} m \ddot{\mathbf{r}} & = m \mathbf{g} - c v^2 \; \hat{\bf{v}} \\ & = m \mathbf{g} - c v \; \mathbf{v} \end{align} $$
The two differential equations for a projectile with quadratic drag
$$ \begin{align} m \dot{v}_x & = -c \sqrt{v_x^2 + v_y^2} \; v_x \\ m \dot{v}_y & = -mg - c \sqrt{v_x^2 + v_y^2} \; v_y \end{align} $$