Letting $g(x) = \sum_0^n a_jx^j$, the coefficient of the $x^{2k}$ term in the product $g(x)g(-x)$ is precisely equal to $$ \alpha_k = \sum_{i+j = 2k} (-1)^j a_i a_j $$ So, if you are given some $f(y) = \sum_{k=0}^n \beta_k y^{k}$, you know that testing whether $f(x^2) = g(x)g(-x)$ reduces to solving the system of multivariate quadratic equations $\alpha_k = \beta_k$ where the $\beta_k$ are specified by your $f$ coefficients and the $\alpha_k$ are as above.
I think no matter which way you cut it, solving a multivariate quadratic system in general is NP hard and unless you get lucky with Buchberger's algorithm, I would not expect an efficient solution.
Update In response to the comment regarding factoring below, yes it is equivalent to testing if $f$ factors. The important thing is: we are testing for factorization rather than computing the factorization. I think the present formulation might offer certain advantages in the case when $f$ does not factor. A certificate of non-existence might be furnished by effective versions of Hilbert's Nullstellensatz.
For instance, if we define the $n+1$ shifted polynomials $$\gamma_k(a_0,\ldots,a_n) = \alpha_k(a_0,\ldots,a_n) - \beta_k,$$ and if there is no set of $a_0,\ldots,a_n$ which simultaneously makes the $\gamma_k$-s vanish, then there must exist polynomials $\delta_k$ so that $$\sum_{k=0}^n \gamma_k\delta_k = 1.$$ More importantly from a computational perspective, the total degree of each $\delta_k$ is bounded by $(n+1)^2 2^{n+1} + 2(n+1)$. For details, see the main theorem of:
Sharp Effective Nullstellensatz. Janos Kollar. Journal of the American Mathematical Society, Vol. 1, No. 4. (Oct., 1988), pp. 963-975