Numerical simulations of evolution of the potentials and impedance spectra of ionselective membranes (ISEs) with ionic sites are presented. The Nernst–Planck–Poisson and continuity equations (NPP) are solved numerically by means of the finite difference method, the Rosenbrock solver and with the use of Matlab platform. Transient solutions for ion-selective electrodes under open- and closed-circuit conditions are computed. The potential-time response to small-current perturbation is used for determination of complex impedances. We present simulations of ISEs as a function of varying diffusivities and ionic concentrations in the “bathing” solutions at interfaces. It is shown that the non-Nernstian behavior of passive membrane electrodes is a result of kinetic constraints at the interfaces, which is manifested in the appearance of an additional arc between the high-frequency bulk and the low-frequency (Warburg) arcs. The presented approach directly relates the diffusivities in the membrane and the interface properties (heterogeneous rate constants determining the transport across interfaces) to the characteristic features of impedance spectra (dimensions and characteristic radial frequencies). NPP problem solved on the Matlab platform allows simulating of the non-linear effects in electro-diffusion.