The simultaneous diffusion of multiple species of dopants in silicon at elevated temperatures is modeled by a nonlinear parabolic system of partial differential equations on a two-dimensional region whose moving boundary depends in part on local dopant concentrations. A numerical solution by finite differences, employing the L-stable TRBDF2 (Trapezoidal Rule with second-order Backward Differentiation Formula) method of time discretization and the "box method" (divergence method) of spatial discretization, is described. Details are given of the methods used to specify curves, to manipulate curves, and to define arbitrary simply-connected regions by their boundary curves. A mathematical analysis of the numerical error of the TRBDF2 method is presented, and numerical experiments comparing alternative choices of matrix representation, timestep size determination, and operator linearization are described and analyzed. |