In this paper we consider systems of weakly coupled nonlinear second order elliptic and parabolic equations with nonlinear, possibly coupled, boundary conditions. The aim is to find invariant sets of the form \[ S = \left\{ {\left( {u_1 ,u_2 , \cdots ,u_m } \right)|\varphi _i (x) \leqq u_i (x) \leqq \psi _i (x){\text{ a.e.}}} \right\}\] for certain nonlinear reaction-diffusion equations, \[ \begin{gathered} U_t + LU = F(U)\quad {\text{in }}\Omega \times [0,\infty ), \hfill \\ BU = G(U)\quad {\text{on }}\partial \Omega \times [0,\infty ), \hfill \\ \end{gathered} \] where $L = (L_1 ,L_2 , \cdots ,L_m )$, ($L_i $ a linear second order elliptic operator), $B = (B_1 ,B_2 , \cdots ,B_m )$, ($B_i $ a linear boundary operator of a general type), and $U = (u_1 ,u_2 , \cdots ,u_m )$. One of the main results says in essence that $S = \{ U\mid \Phi \leqq U \leqq \Psi \} $ is an invariant set if \[L\Phi \leqq F(\Phi )\qquad {\text{and }}L\Psi \geqq F(\Psi )\qquad {\text{in }}\Omega \times [0,\infty ),\] and \[ B\Phi ...