This focus of this paper is the application of the polynomial approximation method to obtain satellite reachability sets for Space Situational Awareness (SSA) applications. Least squares coefficients for the approximation of the final state solution are calculated in a Jacobian free and computationally tractable manner. The numerical integration method used in the polynomial approximation technique exploits the symmetric structure of the input variable probability density functions (pdf’s) to compute the density function for the solution space. The Conjugate Unscented Transformation (CUT) method is utilized for this purpose, and is used to calculate the coefficients with a minimal number of full model propagations. Numerical simulation results are given to validate the approach.