Unbounded solution but bounded Euler discretization Is there an ordinary differential equation in $\mathbb{R}^d$ induced by a gradient vector field with unbounded solutions, for which the difference equations obtained by using the forward Euler method have all solutions bounded?

I read the paper
Guy, Richard T.; Misiurewicz, Michał, Euler approximations can destroy unbounded solutions, Fasc. Math. 44, 43-52 (2010). ZBL1215.65124 (pdf also available here). In this paper, the authors construct examples of differential equations with the form
$$\dot{x}=F(x)$$
that has an unbounded (in forward time) solution $x(t)$, but its Euler approximations for sufficiently small $\epsilon>0$
$$x_n=x_{n−1}+\epsilon F(x_{n−1})$$
have only bounded solutions. However, all examples they propose are not induced by a gradient vector field, i.e., we cannot write $F(x)=\nabla G(x)$ for some $G(x)$, where $\nabla$ stands for the gradient operator. It makes me think whether such phenomenon could happen even for a gradient vector field.

My guess is yes, but it seems hard to find an example that can be verified rigorously by analysis, not just by numerical experiments. Can anyone give a hint on it? It would be better if the example is a polynomial.
 A: The central problem with forward Euler is that it doesn’t continuously update the vector field, and therefore can tunnel through fine features of the vector field.  The gradient case is no exception.  For concreteness, consider the 1D potential $$G(x) =(1/2) x^2 + \sum_{i=1}^{\infty} \int_0^x \eta_{2^{-2-i}}(s - 1 + 2^{-i}) ds$$ where $\eta_{\delta}(\cdot)$ is a smooth nascent delta function with support in $(-\delta/2, \delta/2)$. Note that the supports of these nascent delta functions are disjoint. Starting say from $1/3$, the corresponding continuous solution is infinitely expansive in finite time and hence unbounded; whereas, for any step size $\epsilon>0$, the corresponding Euler solution will tunnel through an infinite # of nascent delta functions and hence remain bounded over a finite time window.
Addendum
In the comments, the OP shared a potential  counterexample, which was brought up in a recent MSE post on August 25th 2022 that remains unanswered. In this "counterexample", the underlying potential/loss/objective/cost function is: $$
G(x,y) = \frac{1}{2} x^2 y^2 - x - y \;.
$$  As noted in the post, based on numerical evidence, it looks like the forward Euler solution is bounded, but the exact solution seems unbounded.  These numerical observations can be validated (or refuted) using a standard stability analysis of the continuous/discrete dynamics outlined below.
First, note that the exact gradient dynamics has only one fixed point at $(1,1)$.  The eigenvalues of the Hessian evaluated at this fixed point $\nabla^2 G(1,1)$ are $3$ and $-1$ with corresponding eigenvectors $(1,1)$ and $(-1,1)$, therefore the fixed point is linearly unstable, and more precisely the instability is in the antidiagonal direction.
To determine the nonlinear stability of the fixed point in the anti diagonal direction, consider the Lyapunov function
$$
\mathcal{V}(x,y) = (x-y)^2 \;,
$$
which measures the distance to the diagonal in the $xy$-plane.  A direct computation shows that
$$
\mathcal{L} V(x,y) = 2 x y \mathcal{V}(x,y) \;,
$$
where $\mathcal{L}$ is the differential operator of the exact gradient dynamics.  This shows that the distance from the diagonal is non-decreasing if $x$ and $y$ have the same sign, and in particular, increasing in a neighborhood of the fixed point $(1,1)$, until either  \begin{align*} 
&y = 1/x^2 \implies \dot{y}=0 ~~\text{and}~~  \dot{x}>0 ~~ \text{or} \\ 
& x=1/y^2 \implies \dot{x}=0 ~~ \text{and} ~~ \dot{y}>0 \;.
\end{align*}
In other words, the exact dynamics becomes (nonlinearly) unstable due to the gravity-like terms in $G(x,y)$, but this instability is weak (non-exponential).
On the other hand, the situation is a bit different along the forward Euler solution.    The fixed point is still unstable with eigenvalues of the linearized (discrete) dynamics being $1-3 \epsilon$ and $1+\epsilon$ with the same eigenvectors.  To assess the nonlinear stability in the anti diagonal direction, we use the same Lyapunov function $$
\mathcal{V}(x_{n+1},y_{n+1}) = (1+\epsilon x_n y_n) \mathcal{V}(x_n, y_n) 
$$
Thus, the distance to the diagonal is increasing in the region
$$
\{ x, y : |1 + \epsilon x y | > 1 \}
$$
which includes the region where $ x y >0$ but also $x y < -2 /\epsilon$.  At a fixed step size, for small time interval the Euler solution essentially behaves like the exact solution, but for larger time intervals instead of asymptoting to one of the following trajectories \begin{align*} 
&y = 1/x^2 \implies \dot{y}=0 ~~\text{and}~~  \dot{x}>0 ~~ \text{or} \\ 
& x=1/y^2 \implies \dot{x}=0 ~~ \text{and} ~~ \dot{y}>0 \;,
\end{align*}
it can catastrophically lose stability and blow up because the component that is small can change sign and the dynamics in that component becomes unstable, and in turn, the dynamics in the other component grows without bound.
This is borne out numerically.  Here is a a trajectory of forward Euler at $\epsilon = 1/6$ in the $xy$-plane over the time interval $[0,10]$ initiated at $(1.6,1.5)$.  The size of the dots increase with time. At first, the Euler solution is essentially following the continuous dynamics and everything looks fine.  However, note that when the $y$-component becomes small, and because forward Euler does not continuously update the gradient, the $y$-component can change sign, which leads to an unexpected/dramatic/artificial instability.

