Let \(\Omega\) be a smoothly bounded domain in \(\mathbb{R}^n\) with unit outward normal \(\eta\), and let \(M_i> 0\), \(\alpha_i\), and \(m_{ij}\) be nonnegative constants for \(i,j=1,\dots, n\). This paper is considered with the initial-boundary value problem \[ \begin{gathered} \frac {\partial u_i}{\partial t} = \nabla(u_i^{\alpha ^i}\nabla u_i) \text{ in } \Omega \times (0,\infty), \\ \frac {\partial u_i}{\partial \eta} = M_i \prod_{j=1}^n u_j ^{m_{ij}} \text{ on } \partial \Omega \times (0,T), \\ u_i(\cdot,0) =u_{i0} \text{ in }\Omega, \end{gathered} \] where \(u_{i0}\) (for \(i=1,\dots,n\)) is a positive \(C^1\) function and \[ \frac {\partial u_{i0}}{\partial \eta} = M_i \prod_{j=1}^n u_{j0} ^{m_{ij}} \text{ on } \partial \Omega \times \{0\}. \] The authors prove a simple existence theorem for this problem. They introduce a nonnegative matrix \(A\), determined explicitly from the constants \(m_{ij}\) and \(\alpha_i\), and show that this problem has a global solution if and only if all the principal minor determinants of \(A\) are nonnegative. The method is based on construction of subsolutions and supersolutions of the system.