I have the following system of PDEs:
$\partial_x q = -c(b_1 \partial_x u + b_2 \partial_y u)$
$\partial_y q = -c(\partial_y u + b_2 \partial_x u)$
where $b_1,b_2,c$ are constants, with $c$ positive. In more compact form:
$\nabla q = -c B \nabla u$
where $B = \left( \begin{matrix} b_1 & b_2 \\ b_2 &1\end{matrix} \right)$.
The matrix $B$ and the function $u$ are given. How can we find $q$?
Lets suppose that $q$ exists. Your right hand side is completely known so you can just integrate. Note that for any $\mathbf x = (x,y)$, $$ \partial_t (q(t\mathbf x)) = \nabla q(t\mathbf x) \cdot \mathbf x$$ so that (assuming some initial condition $q(\mathbf 0)$ is given) $$q(\mathbf x) = q(\mathbf 0) + \int_0^1 \nabla q(t\mathbf x) \cdot \mathbf x \ dt$$ i.e. $$ q(\mathbf x) = q(\mathbf 0) - c\int_0^1 (B\nabla u(t\mathbf x))\cdot \mathbf x \ dt $$
If your domain doesn't contain every straight line from $\mathbf 0$ to any point, a variant of the above will work (since the integral of a gradient does not depend on the path of integration, only the endpoints.)
(edit) OP has clarified that the domain is $\mathbb R^2$. On a simply connected set, a vector field $F$ is a gradient iff it is curl free, so for existence, we will need to assume $$\operatorname{curl} (B \nabla u) = b_2 \partial^2_x u + \partial_{xy} u - b_1 \partial_{xy} u - b_2 \partial_y^2 u = 0.$$ Here, we verify that $\tilde q(\mathbf x) :=q(\mathbf 0) + \int_0^1 F(t\mathbf x)^T \mathbf x dt$, $q(\mathbf 0)$ arbitrary, is indeed a solution to $\nabla q = F$, if $\operatorname{curl}F = \partial_x F_2 - \partial_y F_1 = 0$. Equivalently, the jacobian $\nabla F$ is symmetric.
On the one hand,product and chain rules for differentiating in $\mathbf x$ give $$\partial_i (F(t\mathbf x)_j \mathbf x_j) = \partial_k F_j(t\mathbf x) \partial_i(t\mathbf x)_k\mathbf x_j + F(t\mathbf x)_j \delta_{ij} = t \underbrace{\partial_i F_j(t\mathbf x) \mathbf x_j}_{=(\nabla F(t\mathbf x)^T\mathbf x)_i }+ F_i(t\mathbf x), $$ so that $$ \nabla \tilde q = \int_0^1 t\nabla F(t\mathbf x)^T\mathbf x + F(t\mathbf x) dt.$$
On the other hand, carefully differentiating a specially chosen function in $t$, $$\partial_t(tF(t\mathbf x)_i) = F(t\mathbf x)_i + t\partial_t (F(t\mathbf x)_i) = F(t\mathbf x)_i + t\partial_k F_i(t\mathbf x) \partial_t (t\mathbf x_k) = F(t\mathbf x)_i + t\underbrace{\partial_k F_i(t\mathbf x) x_k}_{=(\nabla F(t\mathbf x) \mathbf x)_i} $$ i.e. $$ \partial_t (tF(t\mathbf x)) = F(t\mathbf x)+t\nabla F(t\mathbf x) \mathbf x$$ which is the integrand of $\nabla \tilde q$, since $\nabla F$ is symmetric; the result now follows from the Fundamental Theorem of Calculus.