Temporal evolution of the energy distributions of fast ions produced by pulses of neutral beam injection (NBI) in fusion devices is studied. Distribution functions of both the NBI ions and the ions born by beam–plasma fusion reactions are considered. New features of these distributions are revealed. In particular, it is found that after a short NBI pulse (blip) the maximum value of the energy distribution function grows or decreases during the slowing down of NBI ions, depending on the particle energy (although the number of particles is conserved); in the first case the distribution looks like a soliton moving in the velocity space. The behaviour is similar for multiple NBI blips that are well separated in time. In contrast to this, when time intervals between the NBI blips are not large, i.e. less than the fast ion slowing down time, a distribution function with several maxima can be formed. In the case of long-lasting NBI pulses, like those in the Wendelstein 7-X stellarator, an inverse energy distribution is formed in a wide region at the initial stage of each pulse, which is a destabilizing factor. The populations of NBI ions and the ions produced by the D–D fusion reaction are evaluated. It is found that the energy content of fusion products is much less than that of the NBI ions for a plasma temperature of about 5 keV and particles injected at 60 keV. Therefore, NBI ions rather than fusion reaction products represent the main source of possible collective phenomena in plasmas under the considered conditions.