In this article, we study square matrices perturbed by a parameter $\epsilon$. An efficient algorithm computing the $\epsilon$-expansion of the eigenvalues in formal Laurent-Puiseux series is provided, for which the computation of the characteristic polynomial is not required. We show how to reduce the initial matrix so that the Lidskii-Edelman-Ma perturbation theory can be applied. We also explain why this approach may simplify the perturbed eigenvector problem. The implementation of the algorithm in the computer algebra system Maple has been used in a quantum mechanics context to diagonalize some perturbed matrices and is available.