We describe an approximate method for treating the mixed quantum-classical (QC) dynamics of many-body systems on N coupled electronic surfaces. The approach is based on calculating NXN reduced Hamiltonian matrices for the classical and quantal degrees of freedom by partial averaging, and then solving the appropriate equations of motion-Hamilton’s equations or the Schrodinger equation-self-consistently. The degrees of freedom requiring a quantum mechanical description are treated using a multistate Schrodinger equation with classically averaged effective time-dependent Hamiltonians and off-diagonal couplings. The classical degrees of freedom are treated by propagating N ensembles of trajectories, one on each electronic surface, using N reduced classical Hamiltonians defined in terms of the expectation value of the full Hamiltonian calculated using the evolving quantum wave functions. An ansatz is proposed to approximately estimate classical off-diagonal density matrix elements required for calculating the classically averaged interactions that couple quantum wave functions on different electronic states. We present the theory and then test it for a simple two-dimensional and two-state model system. Exact quantum and multiconfiguration time-dependent self-consistent-held (MCTDSCF) calculations are carried out to evaluate the QC performance. Good agreement between the MCTDSCF and QC results is obtained for the model considered.