<p>Matrix equations and systems of matrix equations are widely used in control system optimization problems and mathematical economics. However, methods for solving them have only been developed for the most common matrix equations, the Riccati and Lyapunov equations, and there is no universal approach to solving problems of this class. This paper considers methods for solving matrix polynomial equations of arbitrary order with matrix and vector unknowns. An approach to calculating tuples of solutions to polynomial matrix equations based on the theory of branched continued fractions is described. It should be noted that this concerns not only numerical but also symbolic methods of solution. A computational scheme is presented for systems of second-degree polynomial matrix equations with many unknowns. The solution is expanded into a continued matrix fraction. Sufficient conditions for the convergence of continued matrix fractions to solutions and criteria for terminating the calculations in iterative operations are formulated. The results of the numerical experiments confirm the validity of theoretical calculations and the efficiency of the proposed methods.</p>