Electron Cyclotron Resonance (ECR) ion Sources are the most performing machines for the production of intense beams of multi-charged ions in fundamental science, applied physics and industry. Investigation of plasma dynamics in ECRIS still remains a challenge. A better comprehension of electron heating, ionization and diffusion processes, ion confinement and ion beam formation is mandatory in order to increase ECRIS performances both in terms of output beams currents, charge states, beam quality (emittance minimization, beam halos suppression, etc.). Numerical solution of Vlasov equation via kinetic codes coupled to FEM solvers is ongoing at INFN-LNS, based on a PIC strategy. Preliminary results of the modeling will be shown about wave-plasma interaction and electron-ion confinement: the obtained results are very helpful to better understand the influence of the different parameters (especially RF frequency and power) on the ion beam formation mechanism.