This paper presents a modified version of the well-known model of Bernardi and Verbrugge which was developed to simulate the behavior of polymer electrolyte fuel cells. The modified version is based on a different treatment of the electrokinetic model equations, the Butler-Volmer equations. Such equations are analytically integrated in the reactive regions of the electrodes, eliminating the main nonlinear terms in the full mathematical model. It is shown that the modified Bernardi-Verbrugge model is as accurate as the original model, that it allows an extension of the cell current density over which it is possible to find solutions, that the full numerical procedure is very stable, and that the simulations are up to three orders of magnitude faster than those performed with the original model.The mathematical modeling and the numerical simulation of polymer electrolyte membrane fuel cells ͑PEMFCs͒ is receiving, in these last years, more and more attention. Pioneering works can be considered those of Bernardi and Verbrugge 1,2 and those of Springer et al. 3,4 which made a jump from an empirical approach to a more phenomenological and mechanistic type of approach. In these works, the membrane electrode assembly ͑MEA͒ was divided into several regions, and for each region, equations for the charge and mass transport were written down. The description of the MEA was kept at the one-dimensional ͑1-D͒ level, mainly because the essential features of the performance curve of a MEA ͑cell potential vs. cell current density͒ are already well described at this level, while the computational complexity is still affordable.Since then, most of the successive work followed two main streams. In one case, effort has been paid to go from the 1-D modeling to 2-D and even 3-D modeling. These higher dimensionality models offer the possibility of describing the hydrodynamics and multicomponent transport inside the flow channel for reactant distribution ͑see Fig. 1 of the next section͒, and of taking into account different geometries, or flow field configurations, of the flow channels. Clear prototype works of 2-D modeling and simulation have been made by Gurau et al., 5 by Kulikovsky et al., 6 and by Um et al., 7 A 3-D modeling has been developed by Divisek et al., 8 but unfortunately, no implementation of the model has been performed.In the other work stream, attention has been devoted to the inclusion into the models of a more extended phenomenology, such as heat transport, proton transport in membrane with variable water content, membrane permeability, electrode flooding, transient phenomena, etc. ͑for reviews on the state of the art see Gottesfeld and Zawodzinski, 9 and Gurau 10 ͒. However, the approaches developed in these works are, very often, only devoted to partial aspects, losing the ability to construct global models.In this work, we reviewed critically the 1-D model of Bernardi and Verbrugge 1,2 ͑BV model͒, and we developed a modified version of such a model, the modified BV model ͑MBV͒, with the aim to establish a numerical...