In this article we discuss the generalization of a Lagrange multiplier-based fictitious domain method to the simulation of the motion of neutrally buoyant particles in a Newtonian fluid. Then we apply it to study the migration of neutrally buoyant circular cylinders in plane Poiseuille flow of a Newtonian fluid by direct numerical simulation. The Segré–Silberberg effect is found for the cases with one and several circular cylinders. In general, it is believed that the migration away from the center of the channel is due to an effect of the curvature of velocity profile. Via direct numerical simulation, we find that this effect is not weakened by the presence of many particles, but by the collisions among the particles. Experiments and simulations for hundreds of circular cylinder cases show that particles concentrate in the central region where the shear rate is low. A power law associated with the horizontal velocity of the mixture of fluid/particles is also presented.