The present investigation deals with the numerical integration of satellite orbits. As the major contribution to the forces acting on an earth orbiting satellite the anisotropic part of the gravitational field of the earth is given special consideration. The equations of motion are described by a system of first order ordinary differential equations in cartesian earth-fixed coordinates. The main point of the problem is the computation of the gradient of the gravitational potential in the earth-fixed system and the subsequent transformation into the celestial system. The steps which are necessary for this purpose are described in detail. The gravitational potential is approximated by a spherical harmonic expansion of degree and order 360. On this basis a MATLAB program was developed. Several methods for the integration of the equations of motion, available in the MATLAB toolbox, are tested and compared. The comparisons are carried out for two different satellite orbit types, a near-earth orbiting satellite and a high altitude GPS satellite. The final computations are made with the multistep method (predictor-corrector) using the rout ... mehrine ode113 of "The MATLAB Ode Suite". The resulting orbital disturbances are represented as disturbances in the time-dependent corresponding Keplerian orbit elements.