Abstract: The design of three-dimensional borehole trajectories using the natural curve approach under tangent-to-point conditions can be reduced to solving a system of highly nonlinear multivariate equations. Conventional numerical iterative methods often face challenges such as the inability to determine suitable initial iteration values. To address this, a robust and efficient solution strategy is proposed by transforming the multivariate nonlinear system into a univariate nonlinear equation, referred to as the characteristic equation. Once all real roots of the characteristic equation are determined, other variables can be computed through analytical formulas when the inclination angle is prioritized. In cases where the azimuth angle is prioritized, one additional variable is derived by solving a simple trigonometric equation, while the remaining variables are obtained using analytical formulas. The characteristic function is a multimodal continuous function with numerous real roots. To enhance computational efficiency, the constraints are employed to define the maximum permissible interval for meaningful real roots, and the characteristic equation is solved within this interval using root separation and the bisection method. Numerical examples demonstrate that the proposed algorithm effectively solves the borehole trajectory design problem with high speed and robustness, eliminating the need for manually assigned initial iteration values. Moreover, it effectively addresses the issue of multiple solutions in constraint equations. The results are fully consistent with those produced by Compass, a widely used commercial drilling design software. This method has significant practical value and can be integrated into the development of localized drilling design software as an alternative to imported solutions.