The Poisson–Boltzmann (PB) model of electrolyte solutions combines the Poisson equation for electrostatic potentials with a Boltzmann equation c = c0 exp[–eψ/kT] for mobile ion concentrations that is highly nonlinear once the electrostatic potential exceeds 0.1 V. This introduces numerical challenges: first, suitable convergence conditions for the concentration functions become sensitive to the boundary potential. Second, a controlled initial guess must be provided to avoid the finite element method calculation diverging to NaN. We resolve the first challenge by logarithmically scaling the concentration function. A nontrivial log-zero scaling function can handle the near-zero concentrations of coions in a classical point charge model, though is redundant in an advanced model that includes steric forces due to finite ion sizes. The second challenge is resolved with an adaptive throttling algorithm that throttles large values of boundary conditions down to the level of the linear regime and then iteratively raises the throttle until the final nonlinear solution is obtained. The combination of a steric model and throttling enables computation of concentrated electrolytes with electrode potentials as high as 2,000 V. We provide a general derivation of the weak and strong forms of the PB system from the underlying thermodynamic energy functional.