A general mathematical model for the chemical pulping of wood including coupled chemical reactions and diffusion limitations in anisotropic wood chips was developed. The model, which consists of coupled parabolic partial differential equations (PDEs) and ordinary differential equations (ODEs), describes the time-dependent behavior of wood chips as they are exposed to chemicals in the liquid phase. In addition to the reaction−diffusion phenomena, the model describes the change of the wood chip porosity during the process. A numerical algorithm that combines spatial discretization by finite differences with a stiff ODE solver based on the backward difference method was used as an efficient strategy to solve the mass balances of wood chips in batch reactors. Numerical simulations with the software can be used to predict the progress of industrial delignification, that is, production of primarily cellulose through chemical pulping. The effect of the reaction parameters, such as the temperature and the concentrations of the alkaline delignification chemicals, as well as the sizes of the wood chips, can be evaluated with the model, the final goal being the intensification of the chemical pulping process. The model can be used to describe both the current kraft pulping (sulfate pulping) technique, as well as other processes, such as sulfite pulping and pulping in nonaqueous solvents.