(Covers question as originally asked) Well, in the cases you cite, Jordan normal form already takes you quite far. If you want to solve $p(X) = M $ for $n \times n$ complex matrices $X$ and $M$ and $p(t) \in \mathbb{C}[t],$ then if there is a solution at all, the matrix $X$ must act on each generalized eigenspace of $M,$ as it certainly should commute with $M.$ Hence over an algebraically closed field, you may as well reduce to the case that $M$ has a single eigenvalue. Over a field which is not necessarily algebraically closed, you can reduce via rational canonical form to the case that $M$ has a characteristic polynomial which is a power of a singe irreducible polynomial. But it isn't clear to me what generality you want to work in.
Later edit: (Complex case): Solving $p(X) = M$ when $M$ has a single eigenvalue reduces easily to the case when $M$ is nilpotent ( with $p(t)$ replaced by $q(t) - \lambda$ for some scalar $\lambda$). This in turn reduces to piecing together in a consistent fashion solutions to equations $Y^{d} = N^{\prime},$ where both $Y$ and $N^{\prime}$ are nilpotent, and $Y$ has a single Jordan block ( if there is a solution,$M$ must respect any decomposition of the space into indecomposable $X$-invariant summands,and we consider each such indecomposable summand separately). Then it is a question of determining the Jordan normal form of $Y^{d}$ when $Y$ is a nilpotent matrix with a single Jordan block.
Note that if $M$ is diagonalizable, there is always a complex solution, and that, in that case, there are infinitely many solutions if $M$ has an eigenspace of dimension greater than $1$ and $p(t)$ has degree greater than $1.$