Modeling the global distribution of the electric potential in the Earth’s ionosphere is based on the solution of a 2D continuity equation in the ionospheric-magnetospheric current circuit. The potential distribution is described by the boundary value problem for an elliptic system of partial differential equations on the spherical shell approximating the ionosphere, which is divided into three subregions with nonlocal boundary conditions. Implementation of the boundary conditions, which reflect the continuity of the common current circuit and potential equalization at the boundaries of the polar caps, is leading to the mutual dependence of the potential distribution within the northern and southern caps and their influence on the potential distribution in the midlatitude region. The problem is solved by an iterative method with a regularizing operator which is inverted using the separation of the variables and the fast Fourier transform with respect to the azimuthal variable and the sweep method with respect to the latitudinal one.