We present a subdivision based algorithm to compute the solution of an under-constrained piecewise polynomial system of n−2 equations with n unknowns, exploiting properties of B-spline basis functions. The solution of such systems is, typically, a two-manifold in Rn. To guarantee the topology of the approximated solution in each sub-domain, we provide subdivision termination criteria, based on the (known) topology of the univariate solution on the domain’s boundary, and the existence of a one-to-one projection of the unknown solution on a two dimensional plane, in Rn. We assume the equation solving problem is regular, while sub-domains containing points that violate the regularity assumption are detected, bounded, and returned as singular locations of small (subdivision tolerance) size. This work extends (and makes extensive use of) topological guarantee results for systems with zero and one dimensional solution sets. Test results in R3 and R4 are also demonstrated, using error-bounded piecewise linear approximations of the two-manifolds.