In this talk we consider the isothermal Euler-Poisson model with applications to low-temperature low-pressure partially-magnetized plasma discharges. This simple multi-fluid model (ions+electrons) is capable of representing the scale disparity between the different species within plasmas while being theoretically less expensive than kinetic approaches. Nevertheless, the presence of a set of small parameters in the normalized equations, namely the electron-to-ion mass ratio, the spatial scale of electrostatic interactions and the ion-to-electron temperature ratio, present several numerical challenges. It was typically observed that electron current is not accurately predicted in the sheaths (thin boundary layers with charge separation in the vicinity of the walls) by classical approaches, such as second order in time and space Strang operator splitting coupled with MUSCL-Hancock and RK2 methods. Classical method are also constrained, for stability reasons, by the resolution of the small electrostatic spatial scales and their CFL conditions by the characteristic time of evolution of electrons. We propose in this work a series of approaches that can guarantee uniform accuracy across the full domain and/or lift the most restrictive stability conditions. We use the framework of asymptotic-preserving methods in the sense that our methods have convergence properties that remain uniform in all of the asymptotic regimes of the equations (as defined by the set of small parameters) that are present in the considered configurations. We consider implicit-explicit methods, well-balance approaches and high order extensions that allow to gain several order of magnitude in terms of computational cost and render the approach competitive with kinetic models and associated numerical methods. Numerical simulations illustrate the impact of our methods.