In this contribution a numerical model is developed to study the time dynamics of photoluminescence emitted by Tb3+ doped multimode chalcogenide-selenide glass fibers pumped by laser light at approximately 2 microns. The model consists of a set of partial differential equations (PDEs), which describe the temporal and spatial evolution of the photon density and level populations within the fiber. In order to solve numerically the PDEs a Method of Lines is applied. The modeling parameters are extracted from measurements and from data available in the literature. The numerical results obtained support experimental observations. In particular, the developed model reproduces the discrepancies that are observed between the photoluminescence decay curves obtained from different points along the fiber. The numerical analysis is also used to explain the source of these discrepancies.