Find solutions of a dynamical system Let $t\in\mathbb{R}$. Consider the following autonomous dynamical system:
\begin{equation}%\label{eqn: u'(phi)}
x'(t)=-\frac{a_{12}\, r _2-r _1\left(\frac{d_2}{y(t)}+a_{22}\right)}{a_{12}\,a_{21}-\left(\frac{d_1}{x(t)}+a_{11}\right)\left(\frac{d_2}{y(t)}+a_{22}\right)},
\end{equation}
\begin{equation}%\label{eqn: v'(phi)}
y'(t)=-\frac{a_{21}\, r _1-r _2\left(\frac{d_1}{x(t)}+a_{11}\right)}{a_{12}\,a_{21}-\left(\frac{d_1}{x(t)}+a_{11}\right)\left(\frac{d_2}{y(t)}+a_{22}\right)},
\end{equation}
where $a_{ij} (i,j=1,2)$, $d_i (i=1,2)$ and $r_1$ are positive constants; $r_2$ is a negative constant. I have found numerical solutions for some parameters. These numerical solutions show the monotonicity property, i.e. $x'(t)<0$ and $y'(t)>0$. However, I have no idea to give a rigorous proof and determine under which parameters the solutions are monotone. Moreover, I also want to try some standard methods in dynamical systems to find other types of solution which are not monotone. I am not so familiar with dynamical systems. Any suggestion, idea, or comment is welcome, thanks!
 A: First, I propose to rewrite your original system (by multiplying
numerator and denominator by $x\,y$) as
\begin{align}
  x' &= \frac{d_2 r_1\,x + (r_1 a_{22} - r_2 a_{12})\,x\,y}{a_{12}
a_{21}\, x\, y - (d_1 + a_{11}\,x)(d_2 + a_{22}\,y)}, \tag{1a}\\
  y' &= \frac{d_1 r_2\,y + (r_2 a_{11} - r_1 a_{21})\,x\,y}{a_{12}
a_{21}\,x\,y - (d_1 + a_{11}\,x)(d_2 + a_{22}\,y)}. \tag{1b}
\end{align}
If I understand you correctly, you want to see whether the following holds:

CLAIM: Given an arbitrary pair of initial conditions $x(0) =
x_0$, $y(0) = y_0$ such that $x_0 > 0$ and $y_0 < 0$, we have that $x(t)
> 0$ and $y(t) < 0$ for all time $t$.

System $(1)$ is a nonlinear dynamical system, so explicitly solving it
will be difficult. However, we can apply some dynamical systems
techniques to see whether the claim is true.
As a first step, we determine the equilibria of system $(1)$, i.e. the
points $(x,y)$ where $x'=y'=0$. It turns out (I invite you to check
this) that the origin $(x,y)=(0,0)$ is the only equilibrium of the system.
Next, we determine the stability of this equilibrium. We do this by
taking the Jacobian of the right hand side of system $(1)$, and evaluate
it at $(x,y)=(0,0)$. This yields the matrix
\begin{equation}
  J((0,0)) = \begin{pmatrix} r_1/d_1 & 0 \\ 0 & r_2 / d_2
\end{pmatrix}.\tag{2}
\end{equation}
The eigenvalues and eigenfunctions of this matrix can readily be read
off: we have $\lambda_1 = r_1/d_1 > 0$ with eigenvector $(1,0)^T$, and
$\lambda_2 = r_2/d_2 < 0$ with eigenvector $(0,1)^T$. This means (by the
Grobman-Hartman theorem) that we can approximate the phase plane of
system $(1)$ around the origin by the linearised system $(x',y')^T =
J((0,0))(x,y)^T$; from the eigenvalues, we see that the origin is a
saddle. Moreover, the stable manifold of the origin is given by the line
$\{ x=0 \}$, and the unstable manifold of the origin is given by the
line $\{ y=0 \}$. These manifolds act as separatrices in the phase
plane: in other words, because they consist of orbits, other orbits
cannot cross these manifolds and are therefore `caught' between them. In
particular, any orbit which starts in the lower right quadrant will stay
there. So, the claim seems to be true.
However, this linear approximation of the system only holds locally,
that is, sufficiently close to the origin. Because the origin is a
saddle, every orbit (except the ones on the stable manifold) will flow
away from the origin, where the local linear approximation does not
hold anymore. Of course, the stable and unstable manifolds still act as
separatrices, but these are only locally straight. If we zoom out a
little bit, can we still determine what will happen?
Generally, one would now do a so-called manifold expansion, but in the case of system $(1)$, it turns out to be quite easy, because:

OBSERVATION: The line $\{ x=0 \}$ and the line $\{ y=0\}$ are both invariant under the flow of system $(1)$.

So, as you zoom out, you see that the unstable manifold of the origin is exactly equal to the line $\{ y = 0 \}$; also, the stable manifold of the origin is exactly equal to the line $\{x=0\}$. To reiterate, both manifolds consist of orbits, and orbits cannot cross each other, so orbits that are caught in between these two manifolds in the lower right quadrant will stay in the lower right quadrant for all time -- which is exactly the content of the claim.
Obviously, the same statement holds for every other quadrant, so this system indeed preserves some kind of 'monotonicity'.
Question: I'm curious, where does this system come from? It has the form
\begin{equation}
 \vec{x}' = M R \vec{x},
\end{equation}
with $R = \text{diag}(r_1,r_2)$ and $M = \frac{1}{\text{det} B} B^T$, where
$B = S^{-1} \left(D + A\right) S$, with
\begin{equation}
S = \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix} \quad\text{and}\quad A = \begin{pmatrix} a_{11} x & a_{12} x \\ a_{21} y & a_{22} y \end{pmatrix},
\end{equation}
which seems suggestive. 
