Difference equations play a fundamental role in describing the evolution of various phenomena over discrete time intervals. Unlike differential equations, which model continuous processes, difference equations are particularly suitable for systems where changes occur at specific time steps. Over the last twenty-five years of the twentieth century, the theory of difference equations witnessed remarkable development, driven by both theoretical advances and the increasing demand for discrete-time models in applied sciences.

A difference equation is essentially a functional relation among the terms of an unknown sequence, derived from known physical, biological, or economic principles. Once such a relation is established, the equation can be analyzed and solved using various mathematical tools, including analytical methods, numerical approximations, and computational techniques. Of particular interest in modern research are nonlinear and rational difference equations, whose dynamics often exhibit complex behaviors such as oscillations, bifurcations, and chaos.

In this context, several researchers have focused on the qualitative analysis of systems of difference equations. In particular, in 2013, Y. Yazlık, D. T. Tollu, and N. Taşkara investigated the form of solutions of certain systems of rational difference equations. Their work contributed significantly to our understanding of the structural properties and long-term behavior of such systems, and it has motivated further research on stability, boundedness, and convergence properties, especially for systems involving special number sequences such as Fibonacci and Tribonacci numbers.

In this work, we explored the form of solutions, stability character and asymptotics of thenon-linear system of difference equations of the order p + 1.
8><>:
xn+1 = ±1
yn-(p-1)(xn-p±1)+1
x ±1
n-(p-1)(yn-p±1)+1yn+1 = ;
n 2 N0; p ≥ 1
where x
-p; x-(p-1); x-(p-2); : : : ; x0; y-p; y-(p-1); y-(p-2); : : : ; y0 are real initial values with
certain conditions.