**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$.