Recover polynomial \$f(x)\$ from \$f^2(x)\$
Python 2, 638 634 608 603 572 560 508 470 bytes (without brute force as well)
Tthanks to @randomdude999 for the golfing down the code by like 100 bytes as well as everyone who has helped point out smaller golfing tricks!

from gmpy2 import*
R=range
r=reduce
p=lambda n,i=1:n/i and[l+[i]for l in p(n-i,i)]+p(n,i+1)or[[]][n:]
e=lambda f,x:r(lambda y,z:y*x+z,f)
def s(g):
 n=isqrt(len(g));N=n+1;t,_=iroot(abs(g[0]),N)
 for F in(t,-t):
&#9;f=[F]
&#9;for i in R(1,N):f+=[0];f[i]=(g[i]-sum(fac(n)/fac(n-len(j))/r(mul,[fac(j.count(k))for k in R(i+1)])*r(mul,[f[k]for k in j])*F**(N-len(j))for j in p(i))-(i==n)*f[1]*F**~-n)/n/F**n
&#9;if n&lt;2:f[1]=g[1]/-~F
&#9;if all(e(f,e(f,i))==e(g,i)for i in R(N*N)):return f

Try it online!
Explanation
In the code, the functions

m computes the multinomial coefficient \$\binom{n}{k_0,k_1,\dots,n-\sum k}\$
p list out all partitions of n (code from https://codegolf.stackexchange.com/a/84192/61039 (without this code was 601 bytes))
e evaluates our polynoimal
s does the actual functional square root

The way it works is by computing the coefficients of \$f\$ (of degree \$n\$) using the \$n+1\$th highest coefficients of \$f\circ f\$, then checks if the computation is correct as multiple polynomials can lead to the same highest coefficients. The way the computation works is by considering the following series of manipulations:
\begin{align*}
f(x)&amp;=\sum_{i=0}^na_ix^i\\
f(f(x))&amp;=\sum_{i=0}^na_i\left(\sum_{j=0}^na_jx^j\right)^i\\
&amp;=\sum_{i=0}^na_i\sum_{k_0+k_1+\dots+k_n=i}\binom{i}{k_0,k_1,\dots,k_n}\prod_{j=0}^n a_j^{k_j}x^{jk_j}\\
&amp;=\sum_{i=0}^na_i\sum_{k_0+k_1+\dots+k_n=i}\binom{i}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)x^{\sum_{j=0}^njk_j}\\
&amp;=\sum_{k_0+k_1+\dots+k_n\leq n}a_{k_0+k_1+\dots+k_n}\binom{k_0+k_1+\dots+k_n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)x^{\sum_{j=0}^njk_j}\\
&amp;=\sum_{i=0}^{n^2}\left(\sum_{\sum_{j=0}^njk_j=i,k_i\geq0}a_{k_0+k_1+\dots+k_n}\binom{k_0+k_1+\dots+k_n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\right)x^i\\
\end{align*}
From here, if \$i&gt;n^2-n\$, then \$k_1+\dots+k_n=n\$, forcing \$k_0=0\$ since \$a_{&gt;n}=0\$. Furthermore,
$$\sum_{j=0}^njk_j=n\sum_{j=0}^nk_j-\sum_{j=0}^njk_{n-j}=n^2-\sum_{j=0}^njk_{n-j}$$
so the coefficient of \$x^{n^2-i}\$ for \$i&lt;n\$ simplifies to
\begin{align*}
&amp;\hphantom{=}a_n\sum_{\substack{\sum_{j=1}^njk_j=n^2-i\\k_0=0\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\\
&amp;=a_n\sum_{\substack{\sum_{j=0}^njk_{n-j}=i\\\sum_{j=0}^nk_j=n\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_j^{k_j}\right)\\
&amp;=a_n\sum_{\substack{\sum_{j=0}^njk_j=i\\k_0=n-\sum_{j=1}^nk_j\\k_i\geq0}}\binom{n}{k_0,k_1,\dots,k_n}\left(\prod_{j=0}^na_{n-j}^{k_j}\right)
\end{align*}
and the \$\sum_{j=1}^njk_j=i\$ is really asking for partitions of \$i\$.
To extend this analysis to \$i=n\$ to get \$a_0\$, we simply notice that we need to include the case that \$1(n-1)+(n-1)n=n^2\$, i.e. we simply add \$a_0^{n-1}a_1\$.
The reason why this helps is this expansion gives us

\$x^{n^2}\$ coefficient of \$f\circ f\$ is \$a_n^{n+1}\$
\$x^{n^2-i}\$ coefficient is linear in \$a_i\$ for \$0&lt;i&lt;n\$ as the only term containing \$a_i\$ is \$na_0^na_i\$

As all terms in the coefficient of \$x^{n^2-i}\$ in \$f\circ f\$ only depends on \$a_{&lt;i}\$ other than \$na_0^na_i\$, we can easily compute \$a_i\$. The only issue is when \$n\$ is odd, which gives us \$2\$ solutions for \$a_n\$.