We present a generalization of the multiconfigurational time-dependent Hartree (MCTDH) scheme, originally introduced by Meyer, Manthe and Cederbaum [Chem. Phys. Lett. 165, 73 (1990)], to a general nonadiabatic system. In the course of deriving the extended working equations a new compact notation is introduced. Subsequently the equations of motion are applied to a one-dimensional two-surface model system. Calculated energy-resolved transition probabilities for the model system, treated in the MCTDH framework, are shown to be in exact agreement with direct numerically "exact" calculations, using a Split-operator propagation scheme. Finally a comparison is made between the convergence and the consumed CPU-time for the two methods. The two numerical formulations of the scattering problem employ, respectively, a DVR (discrete variable representations) and a FFT (fast Fourier transform) collocation scheme. We also, comment on the use of negative imaginary potentials to remove artificial boundary effects in the two schemes.