  In this paper we study quasilinear elliptic systems given by
  \begin{equation*}
  \begin{aligned}
  -\Delta_{p_1}u_1
  & =-|u_1|^{p_1-2}u_1 \quad && \text{in } \Omega,\newline
  -\Delta_{p_2}u_2
  & =-|u_2|^{p_2-2}u_2 \quad && \text{in } \Omega,\newline
  |\nabla u_1|^{p_1-2}\nabla u_1 \cdot \nu
  &=g_1(x,u_1,u_2) && \text{on } \partial\Omega,\newline
  |\nabla u_2|^{p_2-2}\nabla u_2 \cdot \nu
  &=g_2(x,u_1,u_2) && \text{on } \partial\Omega,
  \end{aligned}
  \end{equation*}
  where $\nu(x)$ is the outer unit normal of $\Omega$ at $x \in
\partial\Omega$, $\Delta_{p_i}$ denotes the $p_i$-Laplacian and $g_i\colon
\partial\Omega \times\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ are
Carath\'{e}odory functions that satisfy general growth and structure conditions
for $i=1,2$. In the first part we prove the existence of a positive minimal and
a negative maximal solution based on an appropriate construction of sub- and
supersolution along with a certain behavior of $g_i$ near zero related to the
first eigenvalue of the $p_i$-Laplacian with Steklov boundary condition. The
second part is related to the existence of a third nontrivial solution by
imposing a variational structure, that is, $(g_1,g_2)=\nabla g$ with a smooth
function $(s_1,s_2)\mapsto g(x,s_1,s_2)$. By using the variational
characterization of the second eigenvalue of the Steklov eigenvalue problem for
the $p_i$-Laplacian together with the properties of the related truncated
energy functionals, which are in general nonsmooth, we show the existence of a
nontrivial solution whose components lie between the components of the positive
minimal and the negative maximal solution.
