In this paper, an advanced algorithm is presented to efficiently form and solve the equations of motion of multibody problems involving uncertainty in the system parameters and/or excitations. Uncertainty is introduced to the system through the application of polynomial chaos expansion (PCE). In this scheme, states of the system, nondeterministic parameters, and constraint loads are projected onto the space of specific orthogonal base functions through modal values. Computational complexity of traditional methods of forming and solving the resulting governing equations drastically grows as a cubic function of the size of the nondeterministic system which is significantly larger than the original deterministic multibody problem. In this paper, the divide-and-conquer algorithm (DCA) will be extended to form and solve the nondeterministic governing equations of motion avoiding the construction of the mass and Jacobian matrices of the entire system. In this strategy, stochastic governing equations of motion of each individual body as well as those associated with kinematic constraints at connecting joints are developed in terms of base functions and modal terms. Then using the divide-and-conquer scheme, the entire system is swept in the assembly and disassembly passes to recursively form and solve nondeterministic equations of motion for modal values of spatial accelerations and constraint loads. In serial and parallel implementations, computational complexity of the method is expected to, respectively, increase as a linear and logarithmic function of the size.