Find $f(f(\cdots f(x)))=p(x)$

$\newcommand{\nest}{\operatorname{nest}}$Let's define a function $\nest(f, x, k)$, which takes a function $f$, an input $x$, and a non-negative integer $k$, and calls $f$ on $x$ repeatedly ($k$ times). For example, $$ \nest(f, x, 0) = x\\ \nest(f, x, 1)=f(x)\\ \nest(f, x, 2)=f(f(x))\\ \nest(f, x, 3)=f(f(f(x)))$$ Formally, this function can be written as $$ \nest(f, x, k)= \begin{cases} x & \text{if } k=0\\ \nest(f, f(x), k-1) & \text{otherwise} \end{cases}$$

For a given $k$ and a polynomial $p$, how can I find a function $f: \mathbb C \to \mathbb C$ such that $\nest(f, x, k)=p(x)$?

If it's not possible to do so in the general case, is it possible with $p(x)=c+x^2$?

You are really asking about functional iterates, normally denoted by $f^ n(x)\equiv \operatorname{nest}(f, x, n)$. I gather you are interested in problem solving methods over rigor.

They are normally obtained by functional iteration through Schröder's equation, but even your simple quadratic $p(x)=x^2+c$ does not have closed solutions, except, e.g., for $c=-2$, a "chaotic" logistic map , as you may see from the examples in the WP article linked above.

In that celebrated special case, a closed form (p302) was found by Ernst Schroeder himself (1870).

Namely, for
$$ p(x)= x^2-2, $$ it follows directly that for $$ y=\frac{x\pm \sqrt{x^2-4}}{2} $$ that is $$ x=y+y^{-1}, $$ one has $$ p(x)=y^2+y^{-2}\equiv p^1 (x). $$ Whence,
$$ p^n (x)= y^{2^n}+ y^{-2^n}. $$

More formally, in E.S.'s language of conjugacy, $$ \psi(p(x))=g(\psi(x)),\\ \psi(x)=\frac{x\pm \sqrt{x^2-4}}{2}\\ g(y)=y^2 \qquad \Longrightarrow \\ g^n(y)=y^{2^n};$$ so that $p(x)= \psi^{-1} \circ g \circ \psi (x)$, and $$p^n= \psi^{-1} \circ g^n \circ \psi ~.$$

I am restricting this to real variables and domains where the objects treated make sense. Your particular question $f^k (x)=x^2-2=p(x)$ then can produce $p^{1/k}(x)$ for suitable domains for you to explore. There are, of course, a plethora of solution-seeking texts on the subject, like C Efthimiou's Introduction to Functional Equations, AMS 2011, ISBN: 978-0-8218-5314-6 , online. A conjugacy iteration approximation method is available in our 2011 paper: Approximate solutions of Functional equations.

I will discuss quadratic $p(x)$, although the general method applies to $p(x)$ with an attractive fixed point. First consider $p(x)=x(a-x)$, with $|a|<1$, then there is an attractive fixed point at $0$. We look for a function $H(x)$ such that equation 1 is obeyed:

$H(p(x))=aH(x)$. ( 1)

We can then put Equation 2:

$p^n(x)=H^{-1} (a^nH(x))$. ( 2)

It is straightforward to expand $H(x)$ as a power series $H(x)= x+x^2/(a(a-1)+...$, but the higher order terms are relatively complex rational functions of a.

One can compute $H(x)$ by computing $y=p^n(x)$ with $n$ sufficiently large that the power series is accurate for $H(y)$ then put $H(x)=a^{-n}H(y)$.

I did numerical experiments for $a=1/2$, and the taylor series for $H(x)$ out to order 10. Equation 1 was obeyed to high accuracy, $H(x)$ appears smooth on its domain, which is the basin of the 0 attractor. $H(x)$ is singular on the boundary of its domain, the Julia set, and I don't think it can be continued beyond there.

$H(x)=0$ for all 0's of $p^n(x)=0$ and any integer $n$. It appears to be analytic on its domain. Equation 2 allows us to differentiate with respect to n and create a continuous flow map in the complex plane. I have not had time to study the inverse map. Clearly it must be multiple valued, but I believe it has regions where single valued definitions exist. In particular, for $p^n(x)=H^{-1} (a^nH(x))$ and $n$ small, we can use $x$ as the starting point of an iterative method to find $H^{-1} (a^nH(x))$.

The function $x^2+c$ has a quadratic fixed point at $\infty$. In this case we transform the mapping to $q(x)=x^2/(1+cx^2)$ and develop a power series about $0$. Because the fixed point is quadratic the functional equation is different. I have not studied this case much yet.