This paper proposes a numerical algorithm designed to integrate the GPS satellite perturbed orbit. The numerical solution is applied to integrate the differential equation of perturbed motion that frames the significant perturbing accelerations. Perturbing potentials are given and the corresponding accelerations in Cartesian coordinates are subsequently deduced. The C++ program that implements a fourth-order Runge-Kutta algorithm is described and comment is made on the perturbed orbit integration results. The paper offers a set of conclusions that will hopefully create a quantitative and qualitative image of GPS orbital perturbations and open a few ways ahead.