The sedimentation of polydisperse suspensions models are defined by systems of N strongly coupled, nonlinear first-order conservation laws, where the solution is a vector of N partial volume fractions or densities. The solution vector should take values in a set of physically relevant values (i.e., the components are nonnegative and sum up at most to a given maximum value). In the 1D case, it is shown that this set, the so-called invariant region, is preserved by numerical solutions produced by a new family of high-order finite volume numerical schemes adapted to this class of models. To achieve this property a pair of linear scaling limiters is applied to a high-order central weighted essentially non-oscillatory (CWENO) polynomial reconstruction to obtain invariant-region-preserving (IRP) high-order polynomial reconstructions. These reconstructions are combined with a first order numerical flux to obtain a high-order numerical scheme for the system of conservation laws. It is proved that this scheme satisfies an IRP property under a suitable CFL condition. For the 2D case, we study a coupled model with a Stokes problem describing the velocity of the mixture and we propose a third order positivity preserving FV-CWENO scheme for the numerical approximation. The theoretical analysis is corroborated with numerical simulations in some scenarios of interest.