A method for high-order treatment of uncertainties in preliminary orbit determination is presented. The observations consist in three couples of topocentric right ascensions and declinations at three observation epochs. The goal of preliminary orbit determination is to compute a trajectory that fits with the observations in two-body dynamics. The uncertainties of the observations are usually mapped to the phase space only when additional observations are available and a least squares fitting problem is set up. A method based on Taylor differential algebra for the analytical treatment of observation uncertainties is implemented. Taylor differential algebra allows for the efficient computation of the arbitrary order Taylor expansion of a suff...