  We present a new computational method for the determination of energy levels
in four-particle systems like H$_2$, HD, and HeH$^+$ using explicitly
correlated exponential basis functions and analytic integration formulas. In
solving the Schr\"odinger equation, no adiabatic separation of the nuclear and
electronic degrees of freedom is introduced. We provide formulas for the
coupling between the rotational and electronic angular momenta, which enable
calculations of arbitrary rotationally excited energy levels. To illustrate the
high numerical efficiency of the method, we present results for various states
of the hydrogen molecule. The relative accuracy to which we determined the
nonrelativistic energy reached the level of $10^{-12}$-$10^{-13}$, which
corresponds to an uncertainty of $10^{-7}$-$10^{-8}$ cm$^{-1}$.
