Consider a real square matrix $A$, lets say $4$ by $4$. Each entry of this matrix $A(i,j)$ is a polynomial function of $x$. For example, $A\_{ij}=a\_{ij}x^3+b\_{ij}x^2+c\_{ij}x+d\_{ij}$. The parameters of $A\_{ij}$, namely $a\_{ij}$,$b\_{ij}$,$c\_{ij}$,$d\_{ij}$, are known explicitly. 


I would like to solve the equation $\det(A(x)) = 0$. It is cumbersome to write down this determinant explicitly as a polynomial in $x$. Is there a way to (1) compute explicitly the coefficients of the scalar polynomial $\det A(x)$, and (2) find its roots?



Numerical ways of solving $\det(A(x)) = 0$
This is called a polynomial eigenvalue problem, and the standard way to solve it is converting it into a standard eigenvalue problem with a technique called [linearization](http://www.maths.manchester.ac.uk/~higham/talks/talk_iwasep06.pdf): if $(Ax^3+Bx^2+Cx+D)v=0$ for some $x\in\mathbb{C}$ and a vector $v\in\mathbb{C}^{n}$, then
$$ \tag{1}
\begin{bmatrix}
B & C & D\\
I \\
& I\\
\end{bmatrix}
\begin{bmatrix}
x^2v\\xv\\v
\end{bmatrix}
=
x
\begin{bmatrix}
A\\&I\\&&I
\end{bmatrix}
\begin{bmatrix}
x^2v\\xv\\v
\end{bmatrix}.
$$
If the matrix $\begin{bmatrix}
A\\&I\\&&I
\end{bmatrix}
$ is nonsingular, then you can multiply both sides by its inverse and reduce to a standard eigenvalue problem; otherwise, it is a problem of the form $Mw=xNw$, for $M,N\in\mathbb{C}^{N\times N},w\in\mathbb{C}^{N}$, which is called generalized eigenvalue problem and can be solved with a [generalization](https://www.mathworks.com/help/matlab/ref/qz.html) of the same [techniques](https://en.wikipedia.org/wiki/QR_algorithm) that are used to solve eigenvalue problems in practice.


This will give you generically $3n$ solutions (or $dn$ for a degree-$d$ polynomial), which are all the solutions of the original problem. Some may appear with higher multiplicity. There are methods to define the geometric multiplicity of an eigenvalue for a matrix polynomial, but they are quite involved (see e.g. the first chapter of *Matrix Polynomials*, Gohberg-Lancaster-Rodman).


If $A$ is singular, there will be fewer than $dn$ solutions; the "missing ones" can be defined to be "eigenvalues at infinity", and this definition behaves consistently with most standard changes of variables and projective geometry tricks (e.g., if you reverse the coefficients of the polynomial you get the inverses of the eigenvalues).


There are alternative ways other than (1) to construct two block matrices to linearize the problem (and determining which one is the best one numerically is still an active research topic); that particular one is known as Frobenius companion form. It is a block version of the well-known [Frobenius companion matrix](https://en.wikipedia.org/wiki/Companion_matrix).


In practice, on a computer, with numerical data, if you have Matlab available, you can just use [`polyeig`](https://it.mathworks.com/help/matlab/ref/polyeig.html?requestedDomain=www.mathworks.com).


If you *really* want the coefficients of the characteristic polynomial, first linearize, then compute a Schur, QZ or Hessenberg factorization of the linearization, and finally compute its characteristic polynomial by evaluation-interpolation. See e.g. [my comment here](https://math.stackexchange.com/questions/405822/what-is-the-fastest-way-to-find-the-characteristic-polynomial-of-a-matrix/405975#comment3382187_405922) for slightly more detail. Note, though, that going through the coefficients of the characteristic polynomial in order to find eigenvalues is not a good idea numerically.


