Numerical solutions to partial differential equations form the backbone of mathematical models that simulate the behavior of various electrochemical systems, specifically, batteries and fuel cells. In this paper, we present a set of numerical algorithms applied to efficiently solve this system of equations. These fast algorithms are identified by fully understanding the physics of the problem and recognizing the strength of the coupling between the governing equations. We illustrate this coupling, specifically in the two potential equations, and demonstrate the need for their simultaneous solution using the Newton method. We take a 2D thermal and electrochemical coupled Li-ion model and extend the familiar Band(J) subroutine by utilizing a Krylov iterative solver, a generalized minimal residual subroutine (GMRES), instead of the direct solver (Gauss elimination), to improve the solution efficiency of the large, nonsymmetric Jacobian system. In addition, we use a nonlinear Gauss-Seidel method to provide the initial guess for the Newton iteration, and precondition the GMRES solver with a block Gauss-Seidel and multigrid algorithm with a smoother based on the tridiagonal matrix algorithm. Every stage in this process has been seen to add to the efficiency of the resulting computer simulation with the final result being a substantial improvement in computation speed, namely, simulating complete discharge of the cell in less than 10 min for grid size of 45 x 32.