In this work, we explore the properties of the matrix method for black hole quasinormal modes on the nonuniform grid. In particular, the method is implemented to be adapted to the Chebyshev grid, aimed at effectively suppressing Runge’s phenomenon. It is found that while such an implementation is favorable from a mathematical point of view, in practice, the increase in precision does not necessarily compensate for the penalty in computational time. On the other hand, the original matrix method, though subject to Runge’s phenomenon, is shown to be reasonably robust and suffices for most applications with a moderate grid number. In terms of computational time and obtained significant figures, we carried out an analysis regarding the trade-off between the two aspects. The implications of the present study are also addressed.