So after reading the derivative refresher document, we come to our
real problem (this is based on heavy digestion of the 
'97 physically-based modeling course notes):

Using the bead on a circular wire example. . .

We have a constraint function of two variables for the
position of the particle. We must take a first and second time
derivative in order to find constraints for velocity and acceleration
of the particle:

C(x,y) = 1/2 * (x^2 + y^2 - 1) = 0

C'(x,y) = 1/2 * (2*x*x' + 2*y*y') = 0
        = x*x' + y*y'

C''(x,y) = (x'*x' + x*x'') + (y'*y' + y*y'') = 0

Legal positions (x,y) must satisfy C(x,y)
Legal velocities (x',y') must satisfy C'(x,y)
Legal accelerations (x'',y'') must satisfy C''(x,y)

If we begin with a legal starting position and velocity
(in other words, if C and C' are satisfied initially), we need
only maintain the acceleration constraint C'' in order to
maintain C and C'. This amounts to adding in a force to correct
for acceleration constraint violations:

  force = mass * acceleration
  F = m * a
  a = F/m

  F = forces acting on particle + constraint correcting force
    = f + f'

  acceleration = (x'',y'') = (f+f')/m = ((fx,fy)+(f'x,f'y))/m
  x'' = (fx+f'x)/m
  y'' = (fy+f'y)/m

Substitute these values back into C'':

  C''(x,y) = (x'*x' + x*(fx+f'x)/m) + (y'*y' + y*(fy+f'y)/m)

We want to find the constraint force f', but we have only
one equation C'' with two unknowns (f'x,f'y). However, since
we require that the constraint force can neither add or remove
energy from the system.

Kinetic energy = 1/2 * mass * Velocity^2
               = 1/2*m*V^2
               = 1/2*m*(x',y')^2
               = 1/2*m*(x'^2 + y'^2)

The time derivative of kinetic energy is the work done =

        = 1/2 * m * (2*x'*x'' + 2*y'*y'')
   work = m * (x'*x'' + y'*y'')

We can substitute the values for x'' and y'' into the work:

        = m * (x'*(fx+f'x)/m + y'*(fy+f'y)/m)
        = x'*(fx+f'x) + y'*(fy+f'y)
        = x'*fx + x'*f'x + y'*fy + y'*f'y
   work = (x'*fx + y'*fy) + (x'*f'x + y'*f'y)

The first term is a function of the forces f acting on the particle.
The second term is a function of the contraint force and should
contribute nothing to the work, so:

  (x'*f'x + y'*f'y) = 0

We will require this equation to be true for every legal velocity.
In other words for every velocity (x',y') that satisfies C':

  (x'*f'x + y'*f'y) = 0
  for all (x',y') such that C'(x,y) = x*x' + y*y' = 0

This simply states that the constraint force f' must point in the
direction of the particle position (x,y). This can be written
by making f' and scalar times the position:

  (f'x,f'y) = lambda * (x,y)

We can substitute this definition of f' into our acceleration
constraint to solve for lambda.

  C''(x,y) = (x'*x' + x*(fx+f'x)/m) + (y'*y' + y*(fy+f'y)/m)
           = (x'*x' + x*(fx+(lambda*x))/m) +
             (y'*y' + y*(fy+(lambda*y))/m)
           = x'*x' + x*fx/m + x^2*lambda/m +
             y'*y' + y*fy/m + y^2*lambda/m
  lambda*(x^2+y^2)/m = -x'*x' - x*fx/m -y'*y' - y*fy/m
  lambda = m * (-x'*x' - x*fx/m - y'*y' - y*fy/m) / (x^2+y^2)
         = (-m*x'^2 - m*y'^2 - x*fx - y*fy) / (x^2+y^2)

We can then use the value for lambda to calculate the constraint
force:

  (f'x,f'y) = lambda * (x,y)

then we can calculate the acceleration:

  x'' = (fx+f'x)/m
  y'' = (fy+f'y)/m

And then we can simply proceed with the normal simulation using
whatever solver we desire.

This calculation is subject to drift we must be corrected by
using a damped spring.