We present a new variational method to study the dynamics of a closed bosonic many-body system, the time-dependent hypernetted-chain Euler-Lagrange method, tHNC. Based on the Jastrow ansatz, it accounts for quantum fluctuations in a non-perturbative way. The tHNC method scales well with the number of dimensions, as demonstrated by our results on one-, two-, and three-dimensional systems. We apply the tHNC method to interaction quenches, i.e. sudden changes of the interaction strength, in homogeneous Bose gases. When the quench is strong enough that the final state has roton excitations (as found and predicted for dipolar and Rydberg-dressed Bose-Einstein condensates, respectively), the pair distribution function exhibits stable oscillations. For validation, we compare tHNC results with time-dependent variational Monte Carlo results in one and two dimensions.