How to find (real-valued) roots of matrix polynomial

Assume you have a fixed ($d=O(1)$ for that matter) degree matrix polynomial $$P(X)=A_0+A_1\cdot X+A_2\cdot X^2+\ldots+A_dX^d$$

Where $A_0,A_1,\ldots A_d\in\mathbb N^{n\times n}$ are given as input. Also given is some constant $\epsilon$.

Can we find a matrix $X_0\in \mathbb R^{n\times n}$ such that $||P(X)||<\epsilon$, or assert such does not exist?

What is the complexity of such algorithm?

There has been a lot of work of numerically finding (approximate) roots of polynomials over the reals, is there an equivalent to matrix polynomials?

For example, if $P(X)=-I+X^2$, then a solution could be $$X= \left( \begin{array}{ccc} 1 & 0 \\ 2015 & -1 \\ \end{array} \right)$$