I have the following system of equations on a 2D domain: \begin{align} \sum_p \mathbf{a}_l^p \cdot \nabla [\alpha g_p] + \beta f_l &= \mathrm{rhs}, \quad l=1, \ldots, N_f,\\ \sum_l \alpha \mathbf{a}_l^p \cdot \nabla f_l + \beta g_p &= 0, \quad p=1, \ldots, N_p, \end{align} where $f_l(\mathbf{x})$ and $g_p(\mathbf{x})$ are the unknowns in some abstract $N_f$- and $N_p$-dimensional spaces ($N_p = N_f + 1$), $\alpha(\mathbf{x})$ and $\beta(\mathbf{x})$ are some functions, and $\mathbf{a}_l^p$ are constant (independent of $\mathbf{x}$, 2D) vectors for each $p$ and $l$ (in particular, $\mathbf{a}_1^1 = \hat{\mathbf{x}}$ and $\mathbf{a}_1^2 = \hat{\mathbf{y}}$).
I want to solve this system using mixed finite elements. After multiplication by test functions $(u, v)$, integrating (by parts in the first equation) and using boundary conditions, I end up with the problem to find $(f, g)$ so that \begin{align} -\sum_p \int \alpha (\mathbf{a}_l^p \cdot \nabla u_l) g_p + \int \beta u_l f_l &= \mathrm{rhs},\\ \sum_l \int \alpha v_p (\mathbf{a}_l^p \cdot \nabla f_l) + \int \beta v_p g_p &= 0 \end{align} for $\forall u, \forall v$.
Note that $\mathbf{a} \cdot \nabla u = \nabla(\mathbf{a} u)$. When $N_f = 1, N_p = 2$, the system reduces to the Poisson-type equation for $f\equiv f_1$ and $\mathbf{g} \equiv (g_1, g_2)$: \begin{align} \nabla(\alpha \mathbf{g}) + \beta f &= \mathrm{rhs},\\ \alpha \nabla f + \beta \mathbf{g} &= 0. \end{align} Here the "charge density" $\mathrm{rhs} - \beta f$ is not constant but depends on $f$. I guess that the well known Raviart-Thomas element should work here. However, I do not see a direct way to generalize this approach to the case $N_f > 1$.
My second intention is to approximate the $H_1$ space. I tried (numerically) the Lagrange elements $P_{k+1}$ for $f_l$ and $P_k$ for $g_p$. For $k=1$ the method fails even for constant $\alpha$ and $\beta$. For $k=2$ and constant $\alpha$ and $\beta$ the solution seems to be quite reasonable, but for some real $\alpha$ and $\beta$ it seems to be unstable.
My question is: how do I choose the correct spaces for the discretization (finite element types)? As far as I understand the method, the spaces should satisfy the LBB conditions, but what the meaningful guesses could be?
Probably, systems of similar kind have already been analyzed in literature. Any references will be highly appreciated.