Summary: We investigate global dynamics of the following systems of difference equations
\[ \begin{cases} x_{n+1} = \frac{\alpha_1 + \beta_1x_n} {A_1+\gamma_n} \\ \gamma_{n+1} = \frac{\gamma_2\gamma_n}{A_2+B_2x_n+\gamma_n}\end{cases} , \quad n=0,1,2,\dots \]
where the parameters \(\alpha_1, \beta_1, A_1, \gamma_2, A_2, B_2\) are positive numbers, and the initial conditions \( x_0\) and \(y_0\) are arbitrary nonnegative numbers. We show that this system has rich dynamics which depends on the region of parametric space. We show that the basins of attractions of different locally asymptotically stable equilibrium points or non-hyperbolic equilibrium points are separated by the global stable manifolds of either saddle points or non-hyperbolic equilibrium points. We give examples of a globally attractive non-hyperbolic equilibrium point and a semi-stable non-hyperbolic equilibrium point. We also give an example of two local attractors with precisely determined basins of attraction. Finally, in some regions of parameters, we give an explicit formula for the global stable manifold.