Multispecies kinematic flow models are defined by systems of strongly coupled, nonlinear first-order conservation laws. They arise in various applications including sedimentation of polydisperse suspensions and multiclass vehicular traffic. Their numerical approximation is a challenge since the eigenvalues and eigenvectors of the corresponding flux Jacobian matrix have no closed algebraic form. It is demonstrated that a recently introduced class of fast first-order finite volume solvers, called polynomial viscosity matrix (PVM) methods [M. J. Castro Diaz and E. Fernandez-Nieto, SIAM J Sci Comput 34 (2012), A2173-A2196], can be adapted to multispecies kinematic flows. PVM methods have the advantage that they only need some information about the eigenvalues of the flux Jacobian, and no spectral decomposition of a Roe matrix is needed. In fact, the so-called interlacing property (of eigenvalues with known velocity functions), which holds for several important multispecies kinematic flow models, provides sufficient information for the implementation of PVM methods. Several variants of PVM methods (differing in polynomial degree and the underlying quadrature formula to approximate the Roe matrix) are compared by numerical experiments. It turns out that PVM methods are competitive in accuracy and efficiency with several existing methods, including the Harten, Lax, and van Leer method and a spectral weighted essentially non-oscillatory scheme that is based on the same interlacing property. (c) 2016 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 32: 1265-1288, 2016