We present an Euler–Lagrange approach for simulating magneto-Archimedes separation of almost neutrally buoyant spherical particles in the flow of a paramagnetic liquid, which is of direct relevance for separating different types of plastic by magnetic density separation. A four-way coupled point-particle method is employed where all relevant interactions between an external magnetic field, a magnetic fluid and discrete immersed particles are taken into account. Particle–particle interaction is modelled by a hard-sphere collision model which takes the interstitial fluid effects into account. First, the motion of rigid spherical particles in a paramagnetic liquid is studied in single- and two-particle systems. We find good agreements between our numerical results and experiments performed in a paramagnetic liquid exposed to a non-homogeneous magnetic field, also in the case of two colliding particles. Next, we investigate the magneto-Archimedes separation of particles with different mass densities in many-particle systems interacting with the fluid. Our results reveal that history effects and interparticle interactions significantly influence the levitation dynamics of particles and have a detrimental impact on the separation performance. We also investigate the effect of particle size and initial distribution on the separation performance. Results show that a reduction in the particle size from 4 to 2 mm leads to a $40\,\%$ increase in the separation time. Moreover, preseparation of particles into two groups of light and heavy particles decreases the separation time by $33\,\%$. The presented method is shown to be a robust and efficient computational framework for the investigation of particle-laden flows of magnetically responsive fluids.