Every day thousands of detections of objects orbiting the Earth are retrieved by observatories. However, a single observation with either radar or optical sensors is in general not sufficient to accurately determine the orbit of an unknown object. Due to the observation geometry and the uncertainty related to sensor accuracy, timing accuracy, and observer state knowledge, we obtain a relatively large orbits set (OS) that is compatible with the single observation rather than a single orbital state. When two independent observations are linked to the same object, the orbit knowledge can be refined by estimating a compatible initial orbit guess from the OSs obtained from the single observations. Our approach for the “linkage problem” is to look at the correlation of the two OSs, which means that they need to be sampled, propagated in time to a common epoch and compared to decide whether they pertain to the same object. The initial effort of our project focuses on the definition of the OS using Differential Algebra (DA). This is obtained by mapping the uncertainties from the measurement space to the state space through a high-order map. This allows us to obtain the object state as a truncated Taylor expansion in the measurements. Once the OSs are defined, the next step is to sample them with the goal of keeping the estimated error introduced by the truncated expansion within a certain threshold in the entire region. The approach for the sampling comes from the estimation of the radius of convergence of the truncated polynomial: outside the convergence disk, the error of the polynomial approximation violates the threshold, which means that most probably the OS cannot be described by just one Taylor expansion but needs to be divided in sub-regions. For this reason, the Automatic Domain Splitting (ADS) tool is introduced. The ADS splits the OS into two or more regions defined by just as many Taylor expansions when the estimated error between the real value and the Taylor approximation exceeds a certain tolerance. This can be done when the OS is created as well as when it is propagated in time, since the error of the Taylor expansion grows along the trajectory due to the nonlinearity of the dynamic. The accuracy of the Taylor expansion is checked by evaluating the error produced by the high-order terms. The region is then split in the direction of the highest estimated error when the threshold is exceeded. The output of the procedure is a set of truncated Taylor expansions that represents the initial OS with the necessary precision in every part of the region. This work presents a comparison between the proposed initial orbit determination approach and alternative methods from the literature to show the improvements we achieved by exploiting accuracy information using DA. The OSs resulting from independent observations are then propagated to a common epoch with the ADS tool and compared to find the probability of intersection, which is used as a measure of the likelihood of their correlation.