Gaussian processes (GP) are becoming more and more popular way to solve statistics and machine learning problems. One of the reasons is the increase in computational power that can handle the inherent computational problem for GP models. Still, in the case of big data, the computational burden can be impractical. For this reason, various approximation methods are developed. In our work, we would like to present an alternative to the internal approximation by using the properties of Chebyshev polynomials. The idea is to calculate the GP model only at Chebyshev nodes and use the property of transforming function values in them to Chebyshev coefficients giving a solution to the original problem. In our research, we propose our version of the algorithm and test it on cases of various functions.