How to solve the following system $\frac{\text{d}x}{\text{d} t} = -Ax + \frac{B}{y} - C$, $ \frac{\text{d}y}{\text{d} t} = -Dx + \frac{E}{y} - F$ Is there a way to analytically solve the following ODE system?
$$
\frac{\text{d}x}{\text{d} t} =
-Ax +  B\left(\frac{1}{y} -1\right) \\
\frac{\text{d}y}{\text{d} t} =
-Cx + D\left(\frac{1}{y} -1\right) 
$$
Where $A,B,C,D>0$ and $x(0)=0,\ y(0)=y_0>0$. $B,D$ may also be treated as including a factor of epsilon, $\varepsilon\ll 1$, however my asymptotics hasn't gone anywhere so far.
 A: There is a constant solution: $y(t) = 1$, $x(t) = 0$.  Moreover, this is a stable equilibrium if $A D - B C > 0$.   
In general you can reduce this system to a single first-order equation for $y$ as a function of $x$.  That differential equation  is, according to Maple, an
Abel's equation of the second kind, class B, and does not seem to have a closed-form general solution.
Of course, in the special case $AD = BC$, you get $\dfrac{dy}{dx} = \dfrac{C}{A}$, and thus $y = y_0 + C x/A$
A: One obtains an Abel differential equation of second kind for $x$ as function of $y$, 
which can be brought into the normal form (changing variables) 
$$
z\frac{dz}{du}-z=f(u),
$$
following standard procedures. In
https://arxiv.org/ftp/arxiv/papers/1503/1503.05929.pdf it is claimed to have solved analitically that normal form (I don't understand fully the article), which would give an exact solution
to the initial problem.
Clearly the equations yield
$$
Axy+By-B=\frac{dx}{dy}(Cxy+Dy-D)
$$
We make the following change of variables: 
$$
z=x+\frac DC-\frac{D}{Cy},
$$
so
$$
x=z-\frac DC+\frac{D}{Cy}\quad\text{and}\quad \frac{dx}{dy}=\frac{dz}{dy}-\frac{D}{Cy^2}.
$$
So the equation reads
$$
Ay\left(z-\frac DC+\frac{D}{Cy}\right)+By-B=\left(\frac{dz}{dy}-\frac{D}{Cy^2}\right)\left(Cy\left(z-\frac DC+\frac{D}{Cy}\right)+Dy-D\right),
$$
and simplifying we obtain
$$
Ayz+\frac{\Delta}{C}(y-1)=\frac{dz}{dy}Cyz-\frac{Dz}{y},\quad\text{where $\Delta=BC-AD$}
$$
which yields
$$
z\frac{dz}{dy}=\left(\frac AC+\frac{D}{Cy^2}\right)z+\frac{\Delta}{C^2}\left(1-\frac 1y\right).
$$
Now we change again the variables, setting $u=\frac AC y-\frac{D}{Cy}$. Then 
$$
\frac{dz}{dy}=\frac{dz}{du}\left(\frac AC+\frac{D}{Cy^2}\right)\quad\text{and}\quad y=\frac{Cu\pm \sqrt{4AD+C^2u^2}}{2A}.
$$
Inserting these values into the equation gives us
$$
\left(\frac AC+\frac{D}{Cy^2}\right)z\frac{dz}{du}=z\frac{dz}{dy}=\left(\frac AC+\frac{D}{Cy^2}\right)z+\frac{\Delta}{C^2}\left(1-\frac 1y\right).
$$
Since 
$$
\frac{\frac{\Delta}{C^2}\left(1-\frac 1y\right)}{\frac AC+\frac{D}{Cy^2}}=\frac{\Delta}{C}\frac{y(y-1)}{Ay^2+D},
$$
we obtain, as announced above,
$$
z\frac{dz}{du}=z+f(u),
$$
where 
$$
f(u)=\frac{\Delta}{C}\frac{y(y-1)}{Ay^2+D}=\frac{\Delta R(R-2A)}{2AC(4AD+CuR)},\quad\text{with $R=Cu\pm \sqrt{4AD+C^2u^2}$}.
$$
A: Depends on what's to be understood by "analytically solve".
One possible approach is via Lie series:

How to properly apply the Lie Series

The main formula in the above reference is:
$$
{\bf x}_1(t) = e^{t X} {\bf x} = e^{t\,{\bf g(x)}\cdot\nabla} {\bf x}
\quad \Longleftrightarrow \quad
\dot{{\bf x}}_1(t) = {\bf g}({\bf x}_1(t)) \quad \mbox{with} \quad {\bf x} = {\bf x}_1(0)
$$
Specified for the OP's case:
$$
\begin{bmatrix} x(t) \\ y(t) \end{bmatrix} = 
e^{t \left\{\left[-Ax_0+B(1/y_0-1)\right]\partial/\partial x_0
+ \left[-Cx_0+D(1/y_0-1)\right]\partial/\partial y_0 \right\} }
\begin{bmatrix} x_0 \\ y_0 \end{bmatrix}
$$
Essentially resulting in a power series expansion:
$$
\begin{bmatrix} x(t) \\ y(t) \end{bmatrix} = \sum_{n=0}^\infty \frac{t^n}{n!}
\left\{\left[-Ax_0+B\left(\frac{1}{y_0}-1\right)\right]\frac{\partial}{\partial x_0}
+ \left[-Cx_0+D\left(\frac{1}{y_0}-1\right)\right]\frac{\partial}{\partial y_0}\right\}^n
\begin{bmatrix} x_0 \\ y_0 \end{bmatrix}
$$
The first few terms are:
$$
\begin{bmatrix} x(t) \\ y(t) \end{bmatrix} = \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} +
t \cdot \begin{bmatrix} -Ax_0+B\left(1/y_0-1\right) \\ -Cx_0+D\left(1/y_0-1\right) \end{bmatrix} + \\
\frac{t^2}{2} \cdot \begin{bmatrix} -(-A x_0 + B (1/y_0 - 1)) A - (-C x_0 + D (1/y_0 - 1)) B / y_0^2 \\
-(-A x_0 + B (1/y_0 - 1)) C - (-C x_0 + D (1/y_0 - 1)) D / y_0^2 \end{bmatrix} + \frac{t^3}{6} \cdots
$$
I'm pretty sure that Robert Israel can provide a little MAPLE program that calculates more of these.
