Modeling of batteries involves simulation of coupled charge and thermal transport, interfacial electrochemical reactions across the porous 3D structure of the electrodes (cathodes and anodes) and electrolyte system. The predominant approach has been the pseudo-2D modeling based on the porous electrode theory where transport in the solid-phase is considered in spherical active particles and the diffusion of Li is solved separately at each point along the thickness direction. These models are limited in their ability to capture spatial variations or able to handle multidimensional structure of the electrodes and do not scale to the system-level. In our approach we adapt the equations based on rigorous homogenization procedures which would allow modeling 3D electrodes. The resulting equations become stiff differential algebraic systems on a single domain and require a robust solution methodology. These formulations are implemented in AMPERES code, an external package for AMP software. In the talk we demonstrate the importance of the coupled solution procedure using Newton methods and accelerate using Krylov iterative solvers. Experiments based on various preconditioning strategies are presented.