The recently introduced full multiple spawning (FMS) method for molecular dynamics beyond the BornOppenheimer approximation is tested against exact numerical solution of the coupled nuclear Schrödinger equation for a two-dimensional model problem with two electronic states. The method uses a multiconfigurational frozen Gaussian ansatz for the wave function and the key idea is to expand the size of the basis set only during nonadiabatic events, using available information to predict the regions of phase space where population will be created. This is accomplished via the spawning procedure which keeps the basis size manageable while ensuring that it provides a reasonable approximation to the exact wave function. The parameters that govern the numerical accuracy of the method are discussed in detail. Expectation values and branching ratios are predicted quantitatively over a broad range of energies.