The Bazant–Storey–Kornyshev (BSK) theory [1–3] recently developed an important continuum framework to expound the nonlocal dielectric permittivity of ionic liquids due to electrostatic correlations, leading to a fourth-order modified Poisson–Fermi equation to model the electrostatic potential field in the solvent (e.g., the electrolyte), while the standard second-order Poisson–Fermi equation is still valid in modeling the electrostatic potential field in the solute. Thus an interface problem is formed between the fourth-order and second-order Poisson–Fermi equations through the interface of solvent and solute with jump coefficients, which has exerted tremendous impacts on applications of electrokinetics, electrochemistry, biophysics, and etc. In this paper, a type of mixed finite element method is developed to solve the proposed interface problem once for all variables: the electrostatic potential, electric field, electric displacement field, electrostatic stress as well as interactional force in the electrolyte in an accurate fashion, and its optimal convergence properties are analyzed for all variables in their respective norms. Numerical experiments are carried out through self-defined mathematical examples to validate all attained theoretical results. Furthermore, as a part of modeling verification for the presented interface problem that models the BSK theory, a practically physical example is investigated to validate the necessity of introducing the fourth-order modified Poisson–Fermi equation to describe the electrostatic correlation effects due to charge reversal phenomenon.