The time period of interest for the dynamic analysis of virtual environments used in haptic rendering often spans one or more changes in contact condition. In this project, we have developed a simulator designed to handle multibody systems with changing constraints, wherein the equations of motion for each mode are formulated in minimal ODE form with constraints embedded before they are passed to an ODE solver. The constraint-embedded equations are formulated symbolically according to a re-combination of terms of the unconstrained equations, and this symbolic process is undertaken on-line by the simulator. Constraint-embedding undertaken on-the-fly enables the simulation of systems with an ODE solver for which constraints are not known prior to simulation start or for which the enumeration of all constraint conditions would be unwieldy because of their complexity or number. Issues of drift associated with DAE solvers that usually require stabilization are sidestepped with the constraint-embedding approach.