  Let $f$ be a polynomial system consisting of $n$ polynomials $f_1,\cdots,
f_n$ in $n$ variables $x_1,\cdots, x_n$, with coefficients in $\mathbb{Q}$ and
let $\langle f\rangle$ be the ideal generated by $f$. Such a polynomial system,
which has as many equations as variables is called a square system. It may be
zero-dimensional, i.e the system of equations $f = 0$ has finitely many complex
solutions, or equivalently the dimension of the quotient algebra $A =
\mathbb{Q}[x]/\langle f\rangle$ is finite. In this case, the companion matrices
of $f$ are defined as the matrices of the endomorphisms of $A$, called
multiplication maps, $x_j : \left\vert \begin{array}{c} h \mapsto x_jh
\end{array} \right.$, written in some basis of $A$. We present a practical and
efficient method to compute the companion matrices of $f$ in the case when the
system is zero-dimensional. When it is not zero-dimensional, then the method
works as well and still produces matrices having properties similar to the
zero-dimensional case. The whole method consists in matrix calculations. An
experiment illustrates the method's effectiveness.
