We report results of systematic simulations of the dynamics of solitons in the framework of the one-dimensional nonlinear Schrödinger equation, which includes the harmonic oscillator potential and a random potential. The equation models experimentally relevant spatially disordered settings in Bose–Einstein condensates (BECs) and nonlinear optics. First, the generation of soliton arrays from a broad initial quasi-uniform state by the modulational instability (MI) is considered following a sudden switch of the nonlinearity from repulsive to attractive. Then, we study oscillations of a single soliton in this setting, which models a recently conducted experiment in a BEC. The basic characteristics of the MI-generated array, such as the number of solitons and their mobility, are reported as functions of the strength and correlation length of the disorder, and of the total norm. For the single oscillating soliton, its survival rate is found. The main features of these dependences are explained qualitatively.