Accurate tracking of satellites that have been and are currently orbitting Earth, has lead to a refined approximation of our planets gravitational potential.
Conversely, the refined planets gravitational potential yield accurate predictions of satellite orbits.

Because of the radial symmetry, the potential function of the gravitational field is formulated as a spherical harmonic expansion with coefficients j, c, and s (see the details below).
The plot aside coutours the potential function at Earths surface level, where rho=r_{E}, and j(2) is set to zero, to enhance details.

We develop a Mathematica notebook, in order to yield the gradient of the potential function algebraically in Euclidean coordinates.
Using the gradient matrix, we employ Runge-Kutta 4-th order integration to compute any satellite orbit in the proximity of Earth.

Remark: MATLAB 7 also supplies a celestial data base - look for "geoid" in the demos.

The gravitational force acting on an object close to Earth, is the negative gradient of the following potential function

evaluated at

The coefficients j, s, c are distilled from long term observations of trajectories of passive satellites that orbit the Earth.
The book Spacecraft System Engineering by Fortescue et al. states the coefficients as

* uses Mathematica's WorldData to plot the countries (instead of a background image), and has fewer example orbits.

Hier stehe ich und kann nicht anders.
Martin Luther

Ground tracks and 3D plots of orbits

We plot several instances of satellite ground tracks.
The points are spaced equidistantly in time and evolve from yellow to red.
Click on the images, for a 3d view of the orbits.

There's nothing out there.
That's why it's called space.
Carl Sagan