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.