In this paper, the Keller box finite-difference scheme is employed in tandem with the so-called boundary immobilization method for the purposes of solving a two-phase Stefan problem that has both phase formation and phase depletion. An important component of the work is the use of variable transformations that must be built into the numerical algorithm in order to resolve the boundary-condition discontinuities that are associated with the onset of phase formation and depletion. In particular, this allows the depletion time to be determined, and the solution to be computed after depletion. The method gives second-order accuracy in both time and space for all variables throughout the entire computation. (C) 2016 Elsevier B.V. All rights reserved.