Heh, it's actually simpler than expected. The system is stable but
not asymptotically stable. Here is a rough sketch of a proof:
First, for simplicity, set
\begin{align*}
\beta &= \frac{b_1}{b_2}\\
\alpha&= \frac{a_1}{a_2}\\
\end{align*}
and then center your system at the equilibrium point $P =
\left(\frac{b_1}{b_2}, \frac{a_1}{a_2}, 0 \right) = \left(\beta,
\alpha, 0 \right)$ which means shift the variables as follows
\begin{align*}
x &\mapsto x + \frac{b_1}{b_2}\\
y &\mapsto y + \frac{a_1}{a_2}\\
z &\mapsto z
\end{align*}
The system becomes
\begin{align*}
\dot{x} &= \beta a_2 \, y + \beta a_3 \, z + a_2
\, xy + a_3 \, xz\\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy\\
\dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz
\end{align*}
which could also be reorganized as
\begin{align*}
\dot{x} &= \beta a_2 \, y + a_2
\, xy + a_3(\beta + x) z\\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy\\
\dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz
\end{align*}
Now observe that the $x,y$-coordinate plane, also given by the
equation $z=0$, is invariant under the flow of the system, meaning
that the vector field is tangent to $z=0$. You can even check
directly that if you take a solution $\big(x(t), y(t) \big)$ of
the 2D system
\begin{align*}
\dot{x} &= \beta a_2 \, y + a_2
\, xy \\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy\\
\end{align*}
then $\big(x(t), \, y(t), \, 0 \big)$ is a solution of the
original system
\begin{align*}
\dot{x} &= \beta a_2 \, y + a_2
\, xy + a_3(\beta + x) z\\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy\\
\dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz
\end{align*}
All of this yields that the plane $z=0$ is a center manifold of
your system and the restriction of the system to that invariant
manifold is
\begin{align*}
\dot{x} &= \beta a_2 \, y + a_2
\, xy \\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy\\
\end{align*} As you already know, the latter is a conservative system
(i.e. integrable, has a first integral) Thus, all trajectories on
$z=0$ are ovals (closed orbits, cycles) and consequently, all
orbits close enough to the invariant plane $z=0$ and in a
neighborhood of the equilibrium point $P$, will tend to approach
it more and more, as long as, the restriction $c_1 <
\frac{b_1}{b_2} c_2$ is satisfied. Thus, they will follow the
cyclic closed orbits on the invariant manifold. Hence, the
equilibrium point is stable, but not asymptotically.
Observe that whenever $c_1 > \frac{b_1}{b_2} c_2$ the equilibrium
$P$ is automatically unstable.
Here is the rigorous proof. The following theorem is a direct corollary from the theory of center manifolds (see for example these slides). One can also prove it, although it requires some techniques from analysis. Notice that in this case, the center manifold is known explicitly and is $z=0$.
Theorem. There exist two small open neighborhood $U, \, V \, \in \, \mathbb{R}^2$ of the point $0 = (0,0,0)$ and a smooth change of variables
$$H \, : \, U \, \to \, V$$
$$H \, : \, (u,v,w) \, \mapsto \, (x,y,z)$$ of the form
\begin{align*}
u &= x + h_1(x,y,z)\\
v &= y + h_2(x,y,z)\\
w &= z
\end{align*}
where the functions $h_1$ and $h_2$ are smooth of order at least two, so that $H$ transforms the system of ODEs
\begin{align}
\dot{u} &= \beta a_2 \, v + a_2
\, uv\\
\dot{v} &= - \alpha b_2\, u - b_2 \, uv \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ (1)\\
\dot{w} &= \left(c_1 - \beta c_2\right) w - g(u,v)w
\end{align}
into the original system of ODEs
\begin{align}
\dot{x} &= \beta a_2 \, y + a_2
\, xy + a_3(\beta + x) z \\
\dot{y} &= - \alpha b_2\, x - b_2 \, xy \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\ (2)\\
\dot{z} &= \left(c_1 - \beta c_2\right) z -c_2 \, xz
\end{align}
Since system (1) is definitely stable, but not asymptotically stable, system (2) is also stable but not asymptotically because (2) is obtained from (1) by the change of variables $H$ and any change of variables preserves stability (which is a topological property).
System (1) is stable because the first two equations are separate from the third. If you pick their first integral, let's call it $J(u,v)$, one can check directly that the function $J(u,v)$ is also a first integral of (1). Then, it is easy to prove the stability of (1). On each level surface $J(u,v) = const$, a general solution $(u(t), v(t), w(t))$ close to $w=0$, which starts from a point on the surface $J(u,v) = const$, stays on the surface $J(u,v) = const$ and because of the coefficient $c_1 - \beta c_2 < 0$ it converges to the solution $(u(t), v(t), 0)$, simply because $\lim_{t \to \infty} \, w(t) = 0$. However,
$(u(t), v(t), 0)$ is actually a closed cycle (a periodic solution). Therefore system (1) is stable, but not asymptotically, because $(u(t), v(t), 0)$ does not converge to the origin.
Best Answer
I give you some hints, I hope it helps. For the lemma: If $\lambda$ is an eigenvalue, then there is a complex eigenvector. Your $u$ and $v$ are exactly the real and imaginary parts of this eigenvectors.
For the Theorem. It is easy if you allow a bit more linear algebra. So assume that everything is complex. Then we can assume without loss of generality that the matrix $A$ is in Jordan normal form.
The solution of the ODE can be represented in general using the matrix exponential function: if your initial value was $x_0\in\mathbb{R}^n$, then the solution wil be $$x(t)=e^{At}x_0.$$ The exponential is usually defined by the power series. Now, if your matrix is in Jordan normal form, then you can calculate this power series. More precisely, if $$A=\begin{pmatrix} \lambda & 1 & 0 & \cdots & 0 \\ 0 & \lambda & 1 & \cdots & 0 \\ \vdots & & & & \vdots \\ 0 &\cdots & &\cdots & \lambda\end{pmatrix} = \text{diag}(\lambda) + N,$$ where $N$ is nilpotent, then $$e^{tN} = \begin{pmatrix} 1 & t & \frac{t^2}{2} & \cdots & \frac{t^{n-1}}{(n-1)!} \\ 0 & 1 & t & \cdots & \frac{t^{n-2}}{(n-2)!} \\ \vdots & & & & \vdots \\ 0 & \cdots & & \cdots & 1 \end{pmatrix}$$, and $$e^{tA} = e^{t\lambda}e^{tN}.$$
You can see now that
I hope this helped.