(Skip the initial two – or even three – paragraphs for the actual solution).

***

The equation $$x_{n+1} - x_n = \nabla F(x_n) - x_n$$ seems to be a discretised version of $$x'(t) = \nabla F(x(t)) - x(t).$$ The vector field $\nabla F(x) - x$ is conservative, with scalar potential $|x|^2/2 - F(x)$ which is confining (that is, it goes to $\infty$ at infinity). Thus, the gradient flow will converge to the set of critical points of the potential. In a typical situation, this set is discrete, in which case the answer for the continuous problem is affirmative (although the limit need not be unique: it may depend on the initial value). If, however, you choose a sufficiently nasty potential, the gradient flow may fail to converge.

***

Now let us try a similar approach for the discrete equation: write $G(x) = |x|^2/2 - \nabla F(x)$, so that $x_{n+1} - x_n = -\nabla G(x_n)$. This is the [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) sequence for the function $G$, with a constant parameter. It is well-known that gradient decent may fail to converge if the parameter is too large, compared to the second derivatives of $G$. However, in our case the second derivative of $G$ in any direction is bounded by $1$, which is precisely what we need.

***

Here the actual solution begins: Write $G(x) = |x|^2/2 - F(x)$, so that $x_{n+1} = x_n - \nabla G(x_n)$. By convexity of $F$,
$$ G(x) \le G(x_n) + \nabla G(x_n) \cdot (x - x_n) + |x - x_n|^2/2 . $$
Setting $x = x_{n+1}$ leads to
$$ \begin{aligned} G(x_{n+1}) & \le G(x_n) + \nabla G(x_n) \cdot (-\nabla G(x_n)) + |\nabla G(x_n)|^2/2 \\ & = G(x_n) - |\nabla G(x_n)|^2/2 . \end{aligned} $$
It follows that $(G(x_n))$ is a (stricly) decreasing sequence, and therefore it converges to a certain limit. In other words, $\sum_n |\nabla G(x_n)|^2 < \infty$, that is, $\sum_{n = 1}^\infty |x_{n+1} - x_n|^2 < \infty$. In particular, $$\lim_{n \to \infty} |\nabla G(x_n)| = \lim_{n \to \infty} |x_{n+1} - x_n| = 0.$$ This implies that the set of accumulation points of $(x_n)$ is connected (see below). Additionally, assuming continuity of $\nabla F$, for every accumulation point $x^* = \lim_{n \to \infty} x_{k_n}$ we have $\nabla G(x^*) = \lim_{n \to \infty} \nabla G(x_{k_n}) = 0$. To summarise:

> If $\nabla F$ is continuous and the set of critical points of $|x|^2/2 - F(x)$ is discrete, then the answer is affirmative: $(x_n)$ converges to a critical point of $F$.

Let me remark why the set of accumulation points of a bounded sequence $(x_n)$, satisfying the condition $\lim_{n \to \infty} |x_{n+1} - x_n| = 0$, is connected. Suppose that it is not, that is, it can be divided into two closed parts $F_1$ and $F_2$, whose distance is positive. Let $3 \delta$ denote the distance between $F_1$ and $F_2$. Choose and arbitrary $N > 0$ large enough, so that $|x_{n + 1} - x_n| < \delta$ when $n > N$. There exist elements $x_i, x_k$ such that $N < i < k$, $\operatorname{dist}(x_i, F_1) < \delta$ and $\operatorname{dist}(x_k, F_2) < \delta$. For $j = i, i+1, \ldots, k$ we have $\operatorname{dist}(x_j, F_1) + \operatorname{dist}(x_j, F_2) \ge 3 \delta$, each summand changes by at most $\delta$ when $k$ is replaced by $k + 1$. It is therefore easy to see that for some $k \in \{i + 1, i + 2, \ldots, k - 1\}$ we must have $\operatorname{dist}(x_k, F_1) \ge \delta$ and $\operatorname{dist}(x_k, F_2) \ge \delta$. In other words, for any $N$ there is $k > N$ such that $\operatorname{dist}(x_k, F_1) \ge \delta$ and $\operatorname{dist}(x_k, F_2) \ge \delta$. By choosing a convergent subsequence of these elements $x_k$, we find an accumulation point of $(x_n)$ away from $F_1 \cup F_2$, a contradiction.

***

OK, I cheated: the actual solution starts here. What happens if the set of critical points of $G(x) = |x|^2/2 - F(x)$ is not discrete? It turns out that $(x_n)$ may fail to converge! It is difficult to describe a counterexample, and I did not check the details; roughly speaking, the graph $G(x)$ described below is an infinite spiral ramp around the unit disk.

We suppose that $d = 2$, we set $c > 0$ small enough (to be determined later), and we identify $(x, y)$ with $x + i y$. The curves $(1 + e^{-t}) e^{i t}$ and $(1 + 2 e^{-t}) e^{i t}$ are [two spirals](http://www.wolframalpha.com/input/?i=ParametricPlot%5B%7B%281%2BE%5E%28-t%29%29%7BCos%5Bt%5D,Sin%5Bt%5D%7D,%281%2B2E%5E%28-t%29%29%7BCos%5Bt%5D,Sin%5Bt%5D%7D%7D,%7Bt,0,4Pi%7D%5D) that will describe the edges of the ramp. The height of the ramp will be roughly $e^{-2 t}$, so the slope will decrease rapidly as the ramp approaches the unit circle.

We define $G(r e^{i t})$ in the following way. First of all, we set $G(r e^{i t}) = 0$ when $r \le 1$. For a fixed $t > 0$, if $1 + e^{-t} \le r < 1 + 2 e^{-t}$, we let $G$ to be a linear function of $r$:
$$ G(r e^{i t}) = c (e^{-2 t} + a(t) (r - 1 - e^{-t})) , $$
where $a(t)$ (extremely small) is chosen in such a way that for every $t > 0$, the point
$$
 G((1 + e^{-t}) e^{i t}) - \nabla G((1 + e^{-t}) e^{i t})
$$
also lies on the same spiral $(1 + e^{-s}) e^{i s}$. A quick-and-dirty calculation suggests that there is a solution $a(t)$ such that the corresponding $s$ satisfies $s \approx t + c e^t a(t)$, and we have $a(t) \sim 2 e^{-3t}$ (and similarly for $a'(t)$ and $a''(t)$); but here I must admit that I did not really check the details.

For $t > 2 \pi$ we define $G(r e^{i t})$ when $1 + 2 e^{-t} < r < 1 + e^{-t + 2 \pi}$ as a cubic polynomial in $r$ that matches the values and the derivatives of $G$ with respect to $r$ at the boundary values $r = 1 + 2 e^{-t}$ and $r = 1 + e^{-t + 2 \pi}$. Another quick-and-dirty calculation shows that wherever we have defined $G$ so far, the second derivatives of $G$ in any direction are bounded by a constant multiple of $c$. We finally extend $G$ to all of $\mathbb{R}^2$ so that $G(r e^{i t}) = 0$ for $r \ge 10$, and the second derivatives of $G$ are everywhere bounded by a constant multiple of $c$.

Now we are ready to eventually choose $c$: we make it small enough, so that the second derivatives of $G$ are in fact bounded by $1$. And we eventually define $F(r e^{i t}) = r^2/2 - G(r e^{i t})$ and $F(r e^{i t}) = F_0(r e^{i t})$ when $r \le 10$.

Now $F$ is convex (because the second derivatives in any direction belong to the interval $(0, 2)$) and $F(r e^{i t}) = r^2 / 2$ for $r \ge 10$. There is one more step needed to assert that $\nabla F$ is bounded: we modify $F(r e^{i t})$ for $r \ge 10$ and we set $F(r e^{i t}) = 10 (r - 10)$ for $r \ge 10$ (this does not affect convexity of $F$).

By the above construction, if $x_0 = 1 + e$, then $x_n = (1 + e^{-t_n}) e^{i t_n}$ for sequence $(t_n)$ such that $t_{n+1} - t_n \to 0$. On the other hand, we must have $t_n \to \infty$, because $\nabla F((1 + e^{-t}) e^{i t}) \ne 0$ for all $t > 0$. Therefore, $(x_n)$ is not convergent.

***

The above construction is not clear, I know. If you are able to simplify it, feel free to edit this answer.