We develop a numerical scheme for Poisson-Boltzmann-Nernst-Planck (PBNP) model. We adopt Gummel's method to treat the nonlinearity of PBNP where Poisson-Boltzmannequation and Nernst-Planckequation are iteratively solved, and then the idea of discontinuous bubble (DB) to solve the Poisson-Boltzmannequation is exploited [6]. First, we regularize the solution of Poisson-Boltzmannequation to remove the singularity. Next, we introduce the DB function as in [6] to treat the nonhomogeneous jump conditions of the regularized solution. Then, we discretize the discontinuous bubble and the bilinear form of Poisson-Boltzmannequation and solve the discretized linear problem by the immersed finite element method. Once Poisson-Boltzmannequation is solved, we apply the control volume method to solve Nernst-Planckequation via an upwinding concept. This process is repeated by updating the previous approximation until the total residual of the system decreases below some tolerance. We provide our numericalexperiments. We observe optimal convergence rates for the concentration variable in all examples having analytic solutions. We observe that our scheme reflects well without oscillations the effect on the distribution of electrons caused by locating the singular charge close to the interface. (C) 2021 Elsevier Inc. All rights reserved.