Given a real function f(X, Y), a box regionB and ε 0, we want to compute an ε-isotopic polygonal approximation to the curve C: f(X, Y) = 0 within B. We focus on subdivision algorithms because of their adaptive complexity. Plantinga & Vegter (2004) gave a numerical subdivision algorithm that is exact when the curve C is non-singular. They used a computational model that relies only on function evaluation and interval arithmetic. We generalize their algorithm to any (possibly non-simply connected) region B that does not contain singularities of C. With this generalization as subroutine, we provide a method to detect isolated algebraic singularities and their branching degree. This appears to be the first complete numerical method to treat implicit algebraic curves with isolated singularities.