This article is concerned with the robust and efficient numerical simulation of case II diffusion, which constitutes an important regime of solvent diffusion into glassy polymers. Even in the one-dimensional case considered here, the numerical simulation of case II diffusion is made difficult by the extreme nonlinearities and coupling in the governing model equations due to a nonlinear flux law needed to produce sharp solvent fronts, a concentration-dependent relaxation time of the polymer used to model the glass-rubber transition, and coupling between the diffusion and deformation phenomena. Having an efficient and accurate solution to such equations is central to advancing a clear understanding of the meaning of such models. The difficulties due to coupling and nonlinearities are highlighted by the consideration of a specific, normalized, one-dimensional case II diffusion model laid out in a general framework of balance laws. Issues such as the stiffness of the spatially discrete differential algebraic equations obtained from the finite element discretization of the governing equations and their bearing on the choice of time-stepping schemes are discussed. The key requirements of numerical schemes, namely, robustness and efficiency, are addressed by the use of an implicit, adaptive, second-order backward differentiation formula with error control for time discretization. Error control is used to maximize the step size to satisfy a target error and the radius of convergence requirements while nonlinear algebraic equations are solved at each time step. An example of an initial boundary value problem is solved numerically to show that the chosen model reproduces case II behavior and to validate that the stated objectives for the numerical simulation are met. Finally, the features and numerical implementation of this model are compared with those of a closely related case II diffusion model due to Wu and Peppas. (C) 2003 Wiley Periodicals, Inc.