An automatic algorithm for accurate numerical gradient calculations has been developed. The algorithm is based on both finite differences and Chebyshev interpolation approximations. The novelty of the method is an automated tuning of the step size perturbation required for both methods. This automation guaranties the best possible solution using these approaches without the requirement of user inputs. The algorithm treats the functions as a black box, which makes it extremely useful when general and complex problems are considered. This is the case of spacecraft trajectory design problems and complex optimization systems. An efficient procedure for the automatic implementation is presented. Several examples based on an Earth-Moon free return trajectory are presented to validate and demonstrate the accuracy of the method. A state transition matrix (STM) procedure is developed as a reference for the validation of the method. / text