This work considers a class of multibody dynamic systems involving bilateral nonholonomic constraints. An appropriate set of equations of motion is employed first. This set is derived by application of Newton's second law and appears as a coupled system of strongly nonlinear second order ordinary differential equations in both the generalized coordinates and the Lagrange multipliers associated to the motion constraints. Next, these equations are manipulated properly and converted to a weak form. Furthermore, the position, velocity and momentum type quantities are subsequently treated as independent. This yields a three-field set of equations of motion, which is then used as a basis for performing a suitable temporal discretization, leading to a complete time integration scheme. In order to test and validate its accuracy and numerical efficiency, this scheme is applied next to challenging mechanical examples, exhibiting rich dynamics. In all cases, the emphasis is put on highlighting the advantages of the new method by direct comparison with existing analytical solutions as well as with results of current state of the art numerical methods. Finally, a comparison is also performed with results available for a benchmark problem.