A new computer program that simulates viscoelastic oxidation of silicon has been developed. Since in this program a tangential procedure is used for time stepping, numerical stability has been improved, and the instability problem that arises from the incorporation of stress dependence into oxide viscosity has been resolved. Thus, oxidation-induced stresses calculated by our program using stress-dependent viscosity have reasonable magnitude over the entire device area. Moreover, in our program, volume expansion due to oxidation of silicon was treated as a dilational strain, as opposed to its treatment as a forced displacement of oxide/silicon interface or a uniaxial strain perpendicular to the interface in most previous programs. Due to this difference, our program can simulate the generation of intrinsic stress below viscous flow temperature (ca. 960°C), and as a result, stress distribution calculated by our program changes drastically at the viscous flow temperature.