Find a formula for the recurrent sequence $q_{n+1}=q_n(q_n+1)+1$ Find an analytic formula for the recurrent sequence $$q_{n+1}=q_n(q_n+1)+1,\;\;q_0\in\mathbb N.$$
(The question was asked on 03.05.2018 by M. Pratsovytyi, see page 109 of Volume 1 of the Lviv Scottish Book). 
 A: If we denote $A_n=q_n+1/2$, then 
$$A_n=A_{n-1}^2+5/4$$
with $A_0=q_0+1/2\ge 3/2$ by $q_0\in\mathbb{N}$.
Further, 
$$\log A_n=2\log A_{n-1}+\log\left(1+\frac{5}{4A_{n-1}^2}\right),$$
namely
$$\frac{1}{2^n}\log A_n-\frac{1}{2^{n-1}}\log A_{n-1}=\frac{1}{2^n}\log\left(1+\frac{5}{4A_{n-1}^2}\right).$$
Thus
$$\log A_n=2^n\left(\log A_{0}+\sum_{k=1}^{n}\frac{1}{2^k}\log\left(1+\frac{5}{4A_{k-1}^2}\right)\right).$$
Clearly,
$$A_n> A_{n-1}^{2^1}> A_{n-2}^{2^2}>\cdots A_0^{2^n}\ge (3/2)^{2^n}.$$
Thus,
$$0<\sum_{k\ge n+1}\frac{2^n}{2^k}\log\left(1+\frac{5}{4A_{k-1}^2}\right)<\sum_{k\ge n+1}\frac{2^n}{2^k}\frac{5}{4A_{k-1}^2}<\frac{5}{4A_{n}^2}.$$
Hence note that $A_0>3/2$ we obtain that
$$1-\frac{5}{4A_n^2}<e^{-\frac{5}{4A_n^2}}\le A_n\kappa^{-2^n}<1$$
with
$$\kappa=A_0\prod_{k=1}^{\infty}\left(1+\frac{5}{4A_{k-1}^2}\right)^{\frac{1}{2^k}}$$ 
a constant depends only on $A_0$. Thus we can prove that
$$\kappa^{2^{n}}-3\kappa^{-2^{n}}<A_n<\kappa^{2^n}$$
for all $n\ge 2$ by note that $\kappa>A_0=3/2$.  For the computing of $\kappa$, it follows from above that
$$A_n^{1/2^n}<\kappa<A_n^{1/2^n}\left(1+\frac{3}{A_n^2}\right)^{1/2^n}<A_n^{1/2^n}\left(1+\frac{3}{2^nA_n^2}\right).$$
A: I'm not quite sure what you mean by an "analytic formula." As Fedor Petrov indicated, there is unlikely to be a closed formula. However, there is a convergent power series. More precisely, consider iteration of the function $f(x)=x^2+x+1$ on the attracting basin at $\infty$, which includes in particular all positive integers. There is an invertible power series, called the Böttcher coordinate of $f$ at $\infty$, that conjugate $f$ to the map $z^2$, and thus it conjugates $f^n$ to the map $z^{2^n}$. This gives an analytic formula is the sense that one gets a convergent power series, but there is not a simple formula for, say, the coefficients of the Böttcher coordinate. Never-the-less, the Böttcher coordinate is a key tool in dynamics for studying iteration in a neighborhood of a superattracting fixed point.
Explicitly, there is an invertible power series
$$\begin{aligned}
\phi(x) &= x + \frac{1}{2} x^2 + \frac{7}{8} x^3 + \frac{3}{4} x^4 + \frac{183}{128} x^5+ O(x^6), \\
\phi^{-1}(x) &= x - \frac{1}{2} x^2 - \frac{3}{8} x^3 + \frac{13}{16} x^4 - \frac{77}{128} x^5+ O(x^6),
\end{aligned}
$$
so that $f(x) = 1/\phi\bigl(\phi^{-1}(1/x)^2\bigr)$. Hence if $q\in\mathbb N$, then
$$
f^n(q) = 1/\phi\bigl(\phi^{-1}(1/q)^{2^n}\bigr).
$$
A: This is a second answer with a somewhat different viewpoint from my other answer. It is an expansion, in some sense, of Gerry Myerson's answer. There is a general theory for estimating these sorts of sequences (even for rational starting points) that goes by the name dynamical canonical heights. I'll restrict attention to this polynomial $f(x)=x^2+x+1$ and integer starting points. For a given starting point $q_0\in\mathbb Z$, the dynamical canonical height of $q_0$ for the map $f$ is given by the limit
$$ \hat h_f(q_0) := \lim_{n\to\infty} \frac{1}{2^n} \log \bigl| f^{\circ n}(q_0)\bigr|.\qquad(*) $$
Here $f^{\circ n}(q_0)$ is the quantity that the OP called $q_n$. It is a standard fact that the limit converges, and that
$$ \log|q_n| = 2^n \cdot \hat h_f(q_0) + O(1), $$
where the $O(1)$ is bounded independently of both $q_0$ and $n$. (It is also relatively easy to give an explicit bound for $O(1)$.) The quantity $\hat h_f(q_0)$ is the logarithm of the $c$ values in Gerry's answer. 
All of this wouldn't be so useful if we had to use the limit formula $(*)$ to compute $\hat h_f(q_0)$, since $(*)$ already requires us to compute $q_n=f^{\circ n}(q_0)$ for large values of $n$. However, there is a rapidly converging series that computes $\hat h_f(q_0)$. It is a modification of a formula originally due to Tate for elliptic curves and can be found in [1]. Computing $k$ terms of this series gives the value of $\hat h_f(q_0)$ with an error of $O(2^{-k})$, so it is quite feasible to compute $\hat h_f(q_0)$ to 100, or even 1000, decimal places.
[1] Call, Gregory S.; Silverman, Joseph H., Canonical heights on varieties with morphisms, Compos. Math. 89, No. 2, 163-205 (1993). ZBL0826.14015.
A: The sequence (at any rate, the case $q_0=1$) has been studied, and references are given at OEIS. The closest thing to a formula given there is $a(n) = [c^{2^n}]$ for $n > 0$, where $c = 1.385089248334672909882206535871311526236739234374149506334120193387331772\dots$
