Asymptotic behavior of system of differential equations Let $f:[0,1]^2 \rightarrow \mathbb{R}^2,$ where $f_1(x,y) = g(y)-x$ and $f_2(x,y) = g(x)-y.$ Here $g(\cdot)$ is a strictly decreasing polynomial function such that $g(0)=1$ and $g(1)=0.$ I am interested in analyzing the asymptotic behavior of the following system of differential equations:
\begin{equation}
\dot{x} = f_1(x,y) \ \text{and} \ \dot{y} = f_2(x,y).
\end{equation}
I want to show that the above system's limiting behavior will be to one of the stationary states i.e., rule out limit cycles and other complicated behavior. Since $\frac{\partial f_1}{\partial x} + \frac{\partial f_2}{\partial y} = -2 \neq 0,$ can I use the Bendixson–Dulac theorem to rule out limit cycles and conclude that the above system converges to one of the stationary states?
The domain in the statement of Bendixson–Dulac's theorem is usually an open set. However my domain of interest $[0,1]^2$ is not open and so am confused. I would really appreciate help on this.
 A: The Bendixson–Dulac theorem is usually stated as you say, but in fact can be proven for arbitrary simply-connected regions (and generalized versions of it also hold for multiply-connected regions). See for instance Theorem 1.8.2 at p. 44 of:
Guckenheimer, J., & Holmes, P. (2013). Nonlinear oscillations, dynamical systems, and bifurcations of vector fields (Vol. 42). Springer Science & Business Media.
Or Problem 7.11 at p. 218 in:
Teschl, G. (2000). Ordinary differential equations and dynamical systems. Graduate Studies in Mathematics, 140, 08854-8019.
A: Preliminary remark. I am not certain whether Bendixon-Dulac grants the global attractiveness of the equilibria. However, via sheer leveraging on the (strict) monotonicity of $g$ and symmetry of the ODE, we can prove that the equilibria is indeed the global attractor which leaves no room for limit cycles or other nontrivial attractors.
The right region in the phase space to be studied is the invariant set $\mathcal{I}=\left\{\left(x,y\right)\in \left[0,1\right]^2\,:\,g(x)\leq y \leq g^{-1}(x)\mbox{ or }g^{-1}(x)\leq y \leq g(x)\right\}$, i.e., the region between the graphs $g$ and $g^{-1}$. This is the right region in the sense that trajectories necessarily accumulate onto it (exponentially fast). Further, observe that $\mathcal{I}$ contains the equilibria given by $\mathcal{E}=\left\{(x,y)\in\left[0,1\right]^{2}\,:\,y=g(x)=g^{-1}(x)\right\}$.
Let $\left(x(t),y(t)\right)_{t\geq 0}$ be the solution to the ODE with initial condition $(x(0),y(0))\in\left[0,1\right]^2$. Define $T$ as the hitting time to hit the set $\mathcal{I}$, i.e., $T\overset{\Delta}=\inf\left\{T\geq 0\,:\, (x(T),y(T))\in \mathcal{I}\right\}$. Let ${\sf d}(w,z)\overset{\Delta}=\|w-z\|_2$.
Claim $1$. If $T=\infty$, then ${\sf d}((x(t),y(t)),\mathcal{E})\overset{t\rightarrow\infty}\longrightarrow 0$ (exponentially fast).

Proof. Define $f_1(t)=\frac{1}{2}\left(x(t)-g(y(t))\right)^2$ and $f_2(t)=\frac{1}{2}\left(x(t)-g^{-1}(y(t))\right)^2$. Assume that $(x(0),y(0))$ is in the lower triangular part of the phase space, i.e., $y(0)<g(x(0))$ and $y(0)<g^{-1}(x(0))$ -- everything that follows applies similarly if we assume an initial condition at the upper triangular part. If $t<T$, then we have that
$$f'_1(t)=(x(t)-g(y(t)))(\dot{x}(t)-g'(y(t))\dot{y}(t))=(x(t)-g(y(t)))(g(y(t))-x(t)-g'(y(t))(g(x(t))-y(t)))\leq -(x(t)-g(y(t)))^2,$$
where the inequality follows since $-g'(y(t))(x(t)-g(y(t)))(g(x(t))-y(t))<0$ as $(g(x(t))-y(t))>0$, $-g'(y(t))>0$ and $(x(t)-g(y(t)))<0$ since $y(t)<g^{-1}(x(t))$.
Therefore, from Grönwall's inequality, $f_1(t)\leq f_1(0) e^{-2t}$. If $T=\infty$, then $f_1(t)\overset{t\rightarrow \infty}\longrightarrow 0$ exponentially fast. Similarly, if $T=\infty$, we conclude that $f_2(t)\overset{t\rightarrow \infty}\longrightarrow 0$ exponentially fast. This is equivalent to ${\sf d}((x(t),y(t)),\mathcal{E})\overset{t\rightarrow \infty}\longrightarrow 0$.

Claim $2$. If $T<\infty$, then ${\sf d}((x(t),y(t)),\mathcal{E})\overset{t\rightarrow\infty}\longrightarrow 0$.

Proof. Now, $(x(T),y(T))\in \mathcal{I}$. Let $g^{-1}(x(T))\leq y(T) \leq g(x(T))$ -- the other case $g(x(T))\leq y(T) \leq g^{-1}(x(T))$ can be dealt with similarly. If $g^{-1}(x(T)) = y(T) = g(x(T))$, then $(x(T),y(T))$ is already at equlibrium. Assume $g^{-1}(x(T))\leq y(T) < g(x(T))$. Remark that $g^{-1}(x(t))\leq y(t) \leq g(x(t))$ for all $t\geq T$. Let us refer to this invariant set as $\mathcal{I}_1$. Further, let $(x^{\star},y^{\star})$ be the equilibrium that lies in the left part of $\mathcal{I}_1$, i.e., $x^{\star}=\sup\limits_{(w_1,w_2)\in \mathcal{E}} w_1 < x(T)$. Then, $V((x,y))\overset{\Delta}= \frac{1}{2}\left(\left(x-x^{\star}\right)^2+\left(y-y^{\star}\right)^2\right)$ conforms to a Lyapunov function granting attractiveness to the equilibrium for any $(x(T),y(T))\in\mathcal{I}_1$: $V$ is definite positive and $\dot{V}(x,y)<0$ for all $(x,y)\in \mathcal{I}_1\setminus \left\{{\sf eq}_{{\sf right}} \right\}$, where ${\sf eq}_{{\sf right}}$ is the equilibrium on the right side. In words, $g^{-1}(x(T))\leq y(T) < g(x(T))$ implies that the solution will acumulate onto the left closest equilibrium, whereas $g(x(T))\leq y(T) < g^{-1}(x(T))$ will imply convergence to the right closest equilibrium.

Claims $1$ and $2$ combined yield the global attractiveness of the equilibria.
Theorem $3$. [$\mathcal{E}$ is the global attractor] ${\sf d}((x(t),y(t)),\mathcal{E})\overset{t\rightarrow\infty}\longrightarrow 0$ regardless of the initial condition $\left(x(0),y(0)\right)\in\left[0,1\right]^2$.
