In fact, the second inequality should be <, giving:
The problem was that on boundaries (when u is an integer), two of the Bi,1(u) functions would evaluate to 1, leading to a total weight of 2 rather than the correct weight of 1.
Once you apply this fix, you will need to do one more thing. The very last iteration on u (where u==umax, where umax is the last value in the knot vector), the point generated will be at (0, 0) because all of the weights are 0. To fix this, do not perform the very last iteration on u, and instead add the last control point directly to the evaluated points.