We study k th order systems of two rational difference equations $$x_n=\frac{\alpha +{\sum }_{i=1}^k{\beta }_i{x}_{n-1}+{\sum }_{i=1}^k{\gamma }_i{y}_{n-1}}{A+{\sum }_{j=1}^k{B}_j{x}_{n-j}+{\sum }_{j=1}^k{C}_j{y}_{n-j}},y_n=\frac{p+{\sum }_{i=1}^k{\delta }_i{x}_{n-i}+{\sum }_{i=1}^k{\epsilon }_i{y}_{n-i}}{q+{\sum }_{j=1}^k{D}_j{x}_{n-j}+{\sum }_{j=1}^k{E}_j{y}_{n-j}}n\in \mathbb{N}$$ . In particular, we assume non-negative parameters and non-negative initial conditions, such that the denominators are nonzero. We develop several approaches which allow us to extend well known boundedness results on the k th order rational difference equation to the setting of systems in certain cases.