  We examine stable solutions of the following symmetric system on a complete,
connected, smooth Riemannian manifold $\mathbb{M}$ without boundary,
\begin{equation*}
  -\Delta_g u_i = H_i(u_1,\cdots,u_m) \ \ \text{on} \ \ \mathbb{M},
  \end{equation*} when $\Delta_g$ stands for the Laplace-Beltrami operator,
$u_i:\mathbb{M}\to \mathbb R$ and $H_i\in C^1(\mathbb R^m) $ for $1\le i\le m$.
This system is called symmetric if the matrix of partial derivatives of all
components of $H$, that is $\mathbb H(u)=(\partial_j H_i(u))_{i,j=1}^m$, is
symmetric. We prove a stability inequality and a Poincar\'{e} type inequality
for stable solutions using the Bochner-Weitzenb\"{o}ck formula. Then, we apply
these inequalities to establish Liouville theorems and flatness of level sets
for stable solutions of the above symmetric system, under certain assumptions
on the manifold and on solutions.
