An algorithm for robust initial orbit determination (IOD) under perturbed orbital dynamics is presented. By leveraging map inversion techniques defined in the algebra of Taylor polynomials, this tool returns a highly accurate solution to the IOD problem and estimates a range centered on the aforementioned solution in which the true orbit should lie. To meet the specified accuracy requirements, automatic domain splitting is used to wrap the IOD routines and ensure that the local truncation error, introduced by a polynomial representation of the state estimate, remains below a predefined threshold. The algorithm is presented for three types of ground-based sensors, namely range radars, Doppler-only radars, and optical telescopes, by considering their different constraints in terms of available measurements and sensor noise. Finally, the improvement in performance with respect to a Keplerian-based IOD solution is demonstrated using large-scale numerical simulations over a subset of tracked objects in low Earth orbit.