The probability and intensity of nuclear reactions involving neutrons are characterized by the corresponding reaction cross-sections which are known to depend strongly on the incident neutron energy. In real applications the neutrons are seldom or never monoenergetic, and are usually characterized by certain continuous energy spectrum. The detailed knowledge of the neutron spectrum is crucial for numerous applications such as the nuclear reactor operation, the traveling wave reactor (TWR) development, including the search of the neutron energy ranges suitable for the wave nuclear burning, the search and prediction of the so-called blowup modes in neutron-multiplying media, the verification of neutron moderation theories and so on. In this paper we describe a method of GEANT4-based Monte Carlo calculation of the neutron spectrum evolution as well as the steady-state neutron spectrum in a system containing a persistent neutron source.