<p>In this study, we investigate the two-dimensional Poisson-Nernst-Planck system using a high-order implicit-explicit backward difference formula for temporal discretization and the integrated radial basis function based on the moving least squares method for spatial discretization on a rectangular domain. This approach results in a sequence of nonlinear algebraic systems. To enhance computational efficiency, we introduce a reduced-order modeling framework by applying proper orthogonal decomposition to the PNP system. However, due to model nonlinearities, the reduced model still exhibits computational complexity dependent on the full system’s variables. To address this, we integrate the discrete empirical interpolation method to further reduce complexity, enabling the expected computational gains from POD while preserving accuracy. Our method captures key physical behaviors of the PNP system, including charge conservation, ionic flux dynamics, and electrostatic interactions. Additionally, it preserves mass conservation and ensures numerical stability, even at high ionic concentrations and strong electric fields. The proposed framework effectively handles nonlinear coupling between diffusion and migration processes, crucial for accurately modeling electrochemical and biological systems. Numerical experiments validate the efficiency and accuracy of the proposed IRBF-MLS-POD-DEIM approach. The results, presented through error analysis, visualizations, and comparative studies, demonstrate that our method not only accelerates computations but also maintains the physical integrity of the PNP system. Additionally, we compare the IRBF-MLS method with its reduced counterparts and the full-order model, highlighting the advantages of the proposed framework in terms of accuracy, execution time, and robustness.</p>