A methodology for initial orbit determination of an object about an astronomical body searches a grid of possible boundary values on the range-to-object over the observation interval to find the grid point and corresponding initial orbit that best fits all of the more than three (N) angles-only observations according to an error metric. The search is conducted by solving a boundary-value problem for different grid points. The state vector is propagated to determine estimated observation directions for the remaining N−2 observations. The grid point (and initial orbit) that best fit all of the observations is selected. For efficiency, the search may start with a broad but coarse initial grid and increase the resolution of the grid as the search converges to a solution.