What is the value of the nested radical $\sqrt[3]{1+2\sqrt[3]{1+3\sqrt[3]{1+4\sqrt[3]{1+\dots}}}}$? The closed-forms of the first three are well-known,
$$x_1=\sqrt{1+\sqrt{1+\sqrt{1+\sqrt{1+\dots}}}}\tag1$$
$$x_2=\sqrt[3]{1+\sqrt[3]{1+\sqrt[3]{1+\sqrt[3]{1+\dots}}}}\tag2$$
$$x_3=\sqrt{1+2\sqrt{1+3\sqrt{1+4\sqrt{1+\dots}}}}\tag3$$
$$x_4=\sqrt[3]{1+2\sqrt[3]{1+3\sqrt[3]{1+4\sqrt[3]{1+\dots}}}}=\;???\tag4$$
with $x_1$ the golden ratio, $x_2$ the plastic constant, and $x_3=3\,$ (by Ramanujan). 

Questions: 



*

*Trying to generalize $x_3$, what is the value of $x_4$ to a $100$ or more decimal places? (The Inverse Symbolic Calculator may then come in handy to figure out its closed-form, if any.)

*What is the Mathematica command to compute $x_4$?


P.S. This other post is related but  only asks for its closed-form which resulted in speculation in the comments. (A method/code to compute $x_4$, and a verifiable numerical value is more desirable.)
 A: We can compute it using backward recursion. Here is a sample Mathematica code: 
nmax = 250; 
x4 = 1; 
Do[ x4 = (1 + i x4)^(1/3) , {i, nmax, 2, -1}]; 
N[x4, 100] 

Mathematica computes this as an expression which it evaluates in the end (i.e. it does not evaluate it to a floating point number in the loop). To guarantee that this matches the true result to $100$ digits one needs to try increasing and increasing nmax untill convergence is seen. For example nmax = 500 and nmax = 250 gives the same $100$ digits.
$$\matrix{
1.7022191326954580969240585907840134288409657961453\\
43207531048888139480023128215942807076912940538302}$$
We can also use this to generate tex-code of the expression. For example taking nmax = 9 and adding x4 // TeXForm gives us this
$$\sqrt[3]{1+2 \sqrt[3]{1+3 \sqrt[3]{1+4 \sqrt[3]{1+5 \sqrt[3]{1+6 \sqrt[3]{1+7
   \sqrt[3]{1+8 \sqrt[3]{10}}}}}}}}$$

Here is some code to evaluate the minimum nmax such that we have convergence to ndigits decimal digits. The idea is to increase nmax by a factor of $2$ until we find convergence and then use bisection on the interval [nmax/2, ..., nmax] to find the minimal value. I'm sure this can be done much better/simpler in Mathematica, but anyway here it is:
(* Evaluates the recursion to level nmax to ndigit precision *)
ExpressionToNdigits[nmax_, ndigits_] := Module[{x4},
   x4 = 1;
   Do[x4 = (1 + i x4)^(1/3), {i, nmax, 2, -1}];
   N[x4, ndigits]
];

(* Finds a nmax such that the recursion have converged to ndigits but it has not converged at nmax/2 *)
FindNmax[ndigits_] := Module[{expression, expressionOld, nmax},
   nmax = 1;
   {expressionOld, expression} = {0, ExpressionToNdigits[nmax, ndigits]};
   While[expression != expressionOld,
    nmax *= 2;
    {expressionOld, expression} = {expression, ExpressionToNdigits[nmax, ndigits]};
   ];
   nmax/2
];

(* Find minimum nmax using bisection *)
FindNmaxMin[ndigits_] := Module[{nmax, nmaxMin, nmaxMax, trueexpression, expression},
   nmaxMax = FindNmax[ndigits];
   nmaxMin = nmaxMax/2;
   trueexpression = ExpressionToNdigits[nmaxMax, ndigits];
   While[nmaxMax - nmaxMin > 1,
    nmax = Floor[(nmaxMin + nmaxMax)/2];
    expression = ExpressionToNdigits[nmax, ndigits];
    If[trueexpression - expression == 0, nmaxMax = nmax, nmaxMin = nmax];
   ];
   nmaxMax
];

FindNmaxMin[100] (* 210 *)
FindNmaxMin[50]  (* 105 *)
FindNmaxMin[20]  (*  42 *)

It turns out the nmax needed to get ndigits precision is almost exactly $n_{\rm max} = 2.1 n_{\rm digits}$.
A: To answer this question we need to find the system that is error resilient.
We can write the equation as
$$y(x)=x\sqrt[3]{1+(x+1)\sqrt[3]{1+(x+2)\sqrt[3]{1+...}}}$$
from where we have
$$y(x)=x\sqrt[3]{1+y(x+1)}$$
or
$$y(x-1)=(x-1)\sqrt[3]{1+y(x)}$$
Now we need to estimate how this function behaves and we can easily see that
$$y(x) \sim x^{\frac{3}{2}}$$
because
$$(x^{\frac{3}{2}})^3\sim(x-1)\sqrt[3]{1+x^{\frac{3}{2}}}$$
From here we have an algorithm. Take large $N$ start with $y(N)=N^{\frac{3}{2}}$ and go backwards using
$$y(k-1)=(k-1)\sqrt[3]{1+y(k)}$$
With, for example, 
$N=20$ we have $y(1)=1.70221913267155$
$N=50$ we have $y(1)=1.70221913269546$ already fixing 14 digits which is obvious 
from $N=100$ $y(1)=1.70221913269546$
It is not difficult to estimate the error terms. If we have missed initial value by $\Delta x$ the error term exponentially diminishes until we reach $y(1)$.
Extracting Mathematica or any other command or language from this is rather trivial.
