Maybe not the answer you are looking for.
To get a more specific answer, more information about the function $f_i$ are
required. In some sense the question is as hard as finding an explicit form for an implicit equation.
Theoretical local state-space form
If you consider the function
$$
g_i(y_i,x_1,x_2) := y_i - f_i(y_i,x_1,x_2)
$$
and let $\widetilde p = (\widetilde y_1, \widetilde x_1, \widetilde x_2)$ be a point with $g(\widetilde p) = 0$.
Then, an application of the implicit function theorem yields a function $\phi_i(x_1,x_2)$ locally defined around $\widetilde p$, such that for all $(y_i,x_1,x_2)$ close to $\widetilde p$ we have
$$
y_i = \phi_i(x_1,x_2) \iff y_i - f_i(y_i,x_1,x_2) = 0.
$$
Therefore, a local state space form is given by $\dot x_i = \phi_i(x_1,x_2)$.
The requirements on $f_i$ are, that $\frac{\partial f_i}{\partial y_i}(\widetilde p) \neq 0$ and that $f_i$ is (at least) continuously differentiable.
Calculus
The derivatives of $\phi_i$ can be expressed in terms of the derivatives of $f_i$,
therefore this approach also allows you to carry out some local analysis.
Theoretical global state-space form
This is in general not possible. Take
$$f(\dot x,x) = (\dot x-1)(\dot x+1),$$
obviously, the local state-space forms are either
$$
\dot x = 1 \quad \text{or} \quad \dot x = -1.
$$
But none of these are global and it depends on the initial conditions which
track you follow. The situation can be arbitrary complex.
Therefore, I would doubt that a general trick exists, which transforms you system
into an explicit equation without strict assumptions on $f_i$.
Numerics
Depending on our task, it might be sufficient to use some root-finding algorithm to get the value of $\phi_i$. (But you need to provide initial guesses for $\dot x_i$, for example the previously computed value of $\phi_i$.)
But, I'm not too familiar with numerics for implicit ODEs.
Best Answer
Hint:
Just notice that: $$(-\dot{f} +af) = -e^{at}(e^{-at}f)'$$ Then, let us denote $$g_1 = e^{-a_1 t}f_1$$ $$g_2 = e^{-a_2 t}f_1$$ The system of two equations is equivalent to $$\begin{align} &\iff \cases{a_1t+\ln(g_1)= A_{11}(a_1t + \ln(-\dot{g_1})) + A_{12}(a_2t + \ln(-\dot{g_2}))\\ a_2t+\ln(g_2)=A_{21}(a_1t + \ln(-\dot{g_1})) + A_{22}(a_2t + \ln(-\dot{g_2}))}\\ &\iff \cases{\ln(g_1)= A_{11} \ln(-\dot{g_1})+ A_{12} \ln(-\dot{g_2})+(A_{11}-1)a_1t + A_{12}a_2t\\ \ln(g_2)= A_{21} \ln(-\dot{g_1}) + A_{22} \ln(-\dot{g_2})+A_{21}a_1t+(A_{22}-1)a_2t }\tag{1} \end{align}$$
We see $(1)$ as a system of linear equations of the two variables $(\ln(-\dot{g_1}),\ln(-\dot{g_2}))$, solve it and we will obtain $$(1)\iff \cases{ \ln(-\dot{g_1})=\alpha_1\ln(g_1)+\beta_1\ln(g_2)+\ln(\gamma_1)t\\ \ln(-\dot{g_2})=\alpha_2\ln(g_1)+\beta_2\ln(g_2)+\ln(\gamma_2)t\\ }\tag{2}$$ where $(\alpha_1,\alpha_2,\beta_1,\beta_2,\gamma_1,\gamma_2)$ are well defined by $(A_{11},A_{12},A_{21},A_{22},a_1,a_2)$
We have: $$(2) \iff \cases{ -\dot{g_1}=\gamma_1^tg_1^{\alpha_1}\cdot g_2^{\beta_1}\\ -\dot{g_2}=\gamma_2^tg_1^{\alpha_2}\cdot g_2^{\beta_2}\\ }\tag{3}$$
Let us define $h_1,h_1$ such that $$ g_1(t)= h_1(\gamma_1^t) \iff h_1(t)= g_1\left(\frac{\ln(t)}{\gamma_1} \right) $$ $$ g_2(t)= h_2(\gamma_2^t) \iff h_2(t)= g_1\left(\frac{\ln(t)}{\gamma_2} \right) $$ then $$(3)\iff \cases{ -\dot{h_1}=\ln(\gamma_1)h_1^{\alpha_1}\cdot h_2^{\beta_1}\\ -\dot{h_2}=\ln(\gamma_2)h_1^{\alpha_2}\cdot h_2^{\beta_2}\\ }\tag{4}$$ we deduce then $$\implies \frac{\dot{h_1}}{\dot{h_2}} = \frac{\ln(\gamma_1)}{\ln(\gamma_2)}h_1^{\alpha_1-\alpha_2}\cdot h_2^{\beta_1-\beta_2}\iff \frac{dh_1}{h_1^{\alpha_1-\alpha_2}}= \frac{\ln(\gamma_1)}{\ln(\gamma_2)} \frac{dh_2}{h_2^{\beta_2-\beta_1}} \tag{5}$$ Integrate $(5)$, we can deduce $$\color{red}{h_2 = \eta \cdot h_1^{\lambda} + \text{const}} \tag{6}$$ Apply $(6)$ to $(4)$, you we have $$\color{red}{-\dot{h_1} = h_1^m(\eta h_1^\lambda + c)^{\beta_2}}$$ and we can determine the function $h_1$ thanks to the hypergeometric function (Ref: WolframAlpha).
After obtaining $h_1$, use $(6)$ we will have $h_2$ and so $(g_1, g_2)$ and $(f_1, f_2)$.