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$?
 A: 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.
A: 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.
A: Similar to tippy2tina's answer, one can easily numerically generate such functions on $\mathbb R$ when they converge monotonically to a fixed-point if a fixed-point exists. Let $f^k(x)=\operatorname{nest}(f,x,k)$. In this case consider for $c\le-1/4$:
$$p(x)=x^2+c,\\q(x)=p^{-1}(x)=\sqrt{x-c}$$
where the inverse, $q$, has $q^n(x)\to\beta=\frac{1+\sqrt{1-4c}}2$ monotonically and
$$q'(\beta)=\frac1{2\beta}$$
describing its rate of convergence
$$q^{n+1}(x)-\beta\sim\frac{q^n(x)-\beta}{2\beta}$$
which implies we can have
$$q^{n+\alpha}(x)-\beta\sim\frac{q^n(x)-\beta}{(2\beta)^\alpha}$$
and thus we may define
$$q^\alpha(x)=\lim_{n\to\infty}p^n\left(\beta+\frac{q^n(x)-\beta}{(2\beta)^\alpha}\right)$$
In your case, you desire to find
$$p^{1/k}(x)=q^{-1/k}(x)=\lim_{n\to\infty}p^n\left(\beta+\frac{q^n(x)-\beta}{(2\beta)^{-1/k}}\right)$$
which will converge to a real solution over $[a,\infty)$ for some $a\le0$. Note that $p$ has two inverses, so only a part of the solution is accounted for, while the other part cannot be real due to issues described in the linked posts.
Graphically for solving in the case of $c=-3.6$ and $k=3:$

