Limit of $(x_n)$ with $0<x_1<1$ and $x_{n + 1} = x_n - x_n^{n + 1}$

Let $0 < x_1 < 1$ and $x_{n + 1} = x_n - x_n^{n + 1}$ for $n \geqslant 1$.

Prove that the limit exists and find the limit in terms of $x_1$.

I have proved the existence but cannot manage the other part.

Thanks for any help.

There is going to be some overlap with @did's post, but what follows is too big to fit into the comment field.

For a fixed, arbitrary $k \geqslant 2$, consider a sequence $y_{n+1} = y_n - y_{n}^{n+k}$ with $y_1 = x_k$. Sequence $\{y_n\}$ is the subsequence of $\{x_n\}$, hence the $x_\ast = \lim_{n \to \infty} x_n$ and $ y_\ast =\lim_{n \to \infty} y_n$ are equal. The limiting value $x_\ast$ is a function of $x_1$, and $y_\ast$ is a function of $x_k$, thus: $$ x_\ast(x_1) = y_\ast(x_k) $$ Few initial terms of the sequence $\{x_n\}$ read: $$ x_2 = x_1 (1-x_1), \quad x_3 = x_1 (1-x_1) \left(1- \left(x_1 (1-x_1)\right)^2 \right) $$

Consider an auxiliary sequence $\{t_{n}\}$ defined as $t_{n+1} = t_n \left( 1- t_1^{n+k} \right)$, $t_1 = y_1$. Because the sequence $y_n$ is decreasing, we have: $$ y_{n+1} = y_n \left( 1- y_n^{n+k-1} \right) \geqslant y_n \left( 1 - y_1^{n+k-1}\right) $$ thus $t_{n} \leqslant y_n$. The sequence $t_n$ is easily solved: $$ t_n = t_1 \prod_{m=1}^{n-1} \left(1-t_1^{m+k-1} \right) = t_1 (t_1^k, t_1)_{n-1} $$ where $(a,q)_n$ denotes the q-Pochhammer symbol.

Therefore, for every $k \geqslant 2$, $$ t_\infty = t_1 \cdot \left(t_1^k, t_1\right)_\infty = x_k \cdot \left(x_k^k, x_k \right)_\infty \leqslant y_\ast = x_\ast $$ The higher the $k$ chosen, the closed $t_\ast$ gets to $x_\ast$.

Here is a numerical verification in Mathematica. Define code for $x_\ast$:

xiter[{xn_, n_}] := {xn - xn^(n + 1), n + 1};
xstar[x1_Real] := 
 First[NestWhile[xiter, {x1, 1}, First[#1] > First[#2] &, 2]]

Now for $y_\ast$:

yiter[k_][{yn_, n_}] := {yn - yn^(n + k), n + 1};
ystar[k_?NumberQ][y1_Real] := 
 First[NestWhile[yiter[k], {y1, 1}, First[#1] > First[#2] &, 2]]

and for $t_\ast$:

tstar[k_?NumberQ][t1_Real] := t1 QPochhammer[t1^k, t1]

For comparison, this is the @did's lower bound:

didstar[x1_Real] := 
 Block[{x2 = x1 (1 - x1)}, x2 (1 - x2 - x2^2)/(1 - x2)]

and code for $x_n$ as a function of $x_1$:

xn[k_Integer][x1_] := 
 First[RecurrenceTable[x[n + 1] == x[n] - x[n]^(n + 1) && x[1] == x1, 
   x, {n, k, k}]]

Here are plots of differences between $x_\ast$ and approximations in natural and logarithmic scales:

You can clearly write each term $x_{n}$ as a polynomial in $x=x_1$, and it should be apparent that all such polynomials are $x + O(x^2)$. Then the difference $x_{n+1}-x_{n}=x_{n}^{n+1}$ is $O(x^{n+1})$, and so the coefficient of $x^{k}$ is the same for all $x_{n}$ with $n\ge k$. This allows us to determine the power series expansion of $x_\infty=\lim_{n\rightarrow\infty}x_{n}$: it is $$ x_\infty(x) = x-x^2-x^3+2x^4+3x^6-20x^7+30x^8-11x^9-31x^{10}+228x^{11}+\dots $$ The OEIS doesn't have anything matching this particular sequence.