A system of rigid bodies with multiple simultaneous unilateral contacts is considered in this paper. The problem is to predict the velocities of the bodies and the frictional forces acting on the simultaneous multicontacts. This paper presents a numerical method based on an extension of an explicit time-stepping scheme and an application of the differential inclusion process introduced by J. J. Moreau. From a differential kinematic analysis of contacts, we derive a set of transfer equations in the velocity-based time-stepping formulation. In applying the Gauss-Seidel iterative scheme, the transfer equations are combined with the Signorini conditions and Coulomb's friction law. The contact forces are properly resolved in each iteration, without resorting to any linearization of the friction cone. The proposed numerical method is illustrated with examples, and its performance is compared with an acceleration-based scheme using linear complementary techniques. Multibody contact systems are broadly involved in many engineering applications. The motivation of this is to solve for the contact forces and body motion for planning the fixture-inserting operation. However, the results of the paper can be generally used in problems involving multibody contacts, such as robotic manipulation, mobile robots, computer graphics and simulation, etc. The paper presents a numerical method based on an extension of an explicit time-stepping scheme, and an application of the differential inclusion process introduced by J. J. Moreau, and compares the numerical results with an acceleration-based scheme with linear complementary techniques. We first describe the mathematical model of contact kinematics of smooth rigid bodies. Then, we present the Gauss-Seidel iterative method for resolving the multiple simultaneous contacts within the time-stepping framework. Finally, numerical examples are given and compared with the previous results of a different approach, which shows that the simulation results of these two methods agree well, and it is also generally more efficient, as it is an explicit method. This paper focuses on the description of the proposed time-stepping and Gauss-Seidel iterations and their numerical implementation, and several theoretical issues are yet to be resolved, like the convergence and uniqueness of the Gauss-Seidel iteration, and the existence and uniqueness of a positive k in solving frictional forces. However, our limited numerical experience has indicated positive answers to these questions. We have always found a single positive root of k and a convergent solution in the Gauss-Seidel iteration for all of our examples.