Husk, 12 bytes

mΣĊ2∂SṪ*z*İ_

Try it online! This uses the standard procedure of calculating -P(x)*P(-x) and dropping the odd-indexed terms. The output coefficients are negated due to the unusual definition of İ_, which is an infinite list of -(-1)^n instead of (-1)^n.

Python + NumPy, 35 bytes

lambda p:p*0+[*p*p(p*0-[1,0])][::2]

Attempt This Online!

Takes and returns instances of numpy.poly1d.

How?

Multplies the given polynomial p(x) with its "mirror" q(x):=p(-x). For each root r of p, q has a root -r and the product has a quadratic factor (x²-r²). Dropping every other coefficient (they are all 0 anyways) in effect substitutes x² by some new variable z, so the resulting polynomial has the correct roots.

The p*0 terms are type coercion starters.

Previous Python + NumPy, 36 bytes

lambda p:p*0+(p*p(p*0-[1,0])).c[::2]

Attempt This Online!

Previous Python + NumPy, 63 bytes

lambda p:P((p*p(-P([1,0]))).c[::2])
from numpy import*
P=poly1d

Attempt This Online!

(Boring) Python + NumPy, 45 bytes

lambda p:poly(roots(p)**2)
from numpy import*

Attempt This Online!

JavaScript (ES7), 60 bytes

This is a port of Neil's answer.

a=>a.map((_,k)=>a.reduce((t,x,i)=>t+(-1)**i*x*~~a[k*2-i],0))

Try it online!

Jelly, 5 bytes

Ær²Æṛ

A monadic Link that accepts and yields coefficients in descending order.

Try it online!

How?

Builtins...

Ær²Æṛ - Link: list of numbers, Coefficients
Ær    - roots of {Coefficients}
  ²   - square {those}
   Æṛ - polynomial with roots {that}

Haskell, 64 bytes

h?(a:b:t)=2*h*b:h?t++[0]
_?l=[0]
f(h:t)=h^2:zipWith(-)(h?t)(f t)

Try it online!

No built-in algebra, just list manipulation. Given \$p\$, we're looking for \$q\$ with

$$ q(x^2) = p(x)p(-x)$$

If we express \$p(x) = h + x p'(x)\$, then we can recurse as

$$p(x)p(-x) = h^2+x^2 \left( 2h \frac{p'(x)-p'(-x)}{2x} - p'(x)p'(-x)\right)$$

So, the first term of \$q(x^2)\$ is \$h^2\$, and the remaining terms are the expression after \$x^2\$. The term \$\frac{p'(x)-p'(-x)}{2x}\$ deletes every other coefficient from \$p'(x)\$. From it, we subtract \$ p'(x)p'(-x)\$, which is the recursive evaluation of our main function on \$p'(x)\$.

If we may stretch polynomial representation to an infinite list trailing with zeroes, we can just write:

Haskell, 51 bytes

h?(a:b:t)=2*h*b:h?t
f(h:t)=h^2:zipWith(-)(h?t)(f t)

Try it online!

Charcoal, 30 24 bytes

ＩＥθΣＥθΣＥθ∧⁼⊗κ⁺μξ×Ｘ±¹ξ×λν

Try it online! Link is to verbose version of code. I/O is a list in ascending order of exponent. Explanation: Directly calculates a polynomial whose roots are the squares of the given polynomial. Works even when the original roots do not have a closed form. Edit: Saved 6 bytes by using @Albert.Lang's optimisation to calculate the sign of each term.

  θ                         Input array
 Ｅ                          Map over values
     θ                      Input array
    Ｅ                       Map over values
        θ                   Input array
       Ｅ                    Map over values
            κ               Outer index
           ⊗                Doubled
          ⁼                 Equals
              μ             Inner index
             ⁺              Plus
               ξ            Innermost index
         ∧                  Logical And
                   ¹        Literal integer `1`
                  ±         Negated
                 Ｘ          Raised to power
                    ξ       Innermost index
                ×           Multiplied by
                      λ     Inner value
                     ×      Multiplied by
                       ν    Innermost value
      Σ                     Take the sum
   Σ                        Take the sum
Ｉ                           Cast to string
                            Implicitly print

Wolfram Language (Mathematica), 21 bytes

Resultant[#,x^2-a,x]&

Try it online!