Density functional theory (DFT) offers computationally affordable way of describing static and dynamic properties of superfluid ⁴He. In general, the DFT models yield single particle-like Schrödinger equations with a nonlinear potential term that accounts for all the many-body interactions. The resulting equations can be solved for small amplitude plane wave excitations in the bulk whereas fully numerical solution must be sought in more complicated cases. In this paper we propose a numerical method that can be used in solving the time-dependent nonlinear Schrödinger equation in both real and imaginary times. The method is based on operator splitting technique where each component operator is treated with a unitary semi-implicit Crank–Nicolson scheme. In order to increase the stability of the method for complex valued nonlinear potentials, a predict–correct scheme is employed in the simulations. The numerical calculations indicate that the scheme is numerically sufficiently stable and well behaving, exhibits high degree of parallelism, and produces results in agreement with the existing numerical data. In the numerical examples we apply the method to obtain dispersion relation of the bulk liquid and to calculate solvation and absorption spectrum of atomic boron solvated in superfluid helium.