Numerical simulation of heterojunction solar cells is a useful approach to investigate the inner physics of the devices, and can be implemented by simulation tools today. An intra-band tunneling model has been described previously and applied to some heterojunction devices. However, the commonly used algorithm for solving the associated equations, the Gummel method, fails in the case of high recombination rates such as in amorphous Si solar cells. In this work, we present an improved algorithm that enhances convermicmgence when solving for intra-band tunneling in heterojunction solar cells with high defect state densities. The algorithm uses a sequential application of the Gummel and Newton methods tailored for optimal stability and convergence. As an example, simulation results for a heterojunction with intrinsic thin layer solar cell are discussed. The improved algorithm and related simulations in this work are implemented through a program we have developed, wxAMPS. Crown Copyright (c) 2012 Published by Elsevier B.V. All rights reserved.