The vortex dynamics in mesoscopic superconducting cylinders with rectangular cross section under an axially applied magnetic field is investigated in the multivortex London regime. The rectangles considered range from a square up to an infinite slab. The flux distribution and total flux carried by a vortex placed in an arbitrary position of the sample is calculated analytically by assuming Clem's solution for the vortex core. The Bean-Livingston energy barrier is also analytically calculated in this framework. A Langevin algorithm simulates the flux penetration and dynamical evolution of the vortices as the external field is slowly cycled. The simulated magnetization process is governed by metastable states. The magnetization curves are hysteretic, with paramagnetic response in part of the downward branch, and present a series of peaks corresponding to the entry or expulsion of a single vortex. For elongated rectangles, the vortices arrange themselves into parallel vortex chains and an additional modulation of the magnetization, corresponding to creation or destruction of a vortex chain, comes out.