New numerical procedures are proposed to solve the symmetric matrix polynomial equation A(T)(-s) X(s) + X-T(-s) A(s)= 2B(s) that is frequently encountered in control and signal processing. An interpolation approach is presented that takes full advantage of symmetry properties and leads to an equivalent reduced-size linear system of equations. It results in a simple and general characterization of all solutions of expected column degrees. Several new theoretical results concerning stability theory and reduced Sylvester resultant matrices are also developed and used to conclude a priori on the existence of a solution. By means of numerical experiments, it is shown that our algorithms are more efficient than older methods and, namely, appear to be numerically reliable.