This paper presents a numerical method for deriving a symplectic state transition matrix for high-fidelity Earth orbits subject to non-dissipative perturbation forces. By taking advantage of properties of Hamiltonian systems, this method provides an exact solution space mapping of linearized orbital dynamics, preserving the symplectic structure that all Hamiltonian systems should possess by nature. This method can be applied to accurate, yet computationally efficient dynamic filters, long-term propagations of the motions of formation flying spacecraft and the eigenstructure analysis of N-body dynamics, etc., when the exact structure-preserving property is crucial. We show the derivation of the numerical method of symplectic state transition matrix, and apply it to Earth orbits with perturbation forces based on real ephemerides. These numerical examples reveal that this method shows improvements in preserving the structural properties of the state transition matrix, and in the computational efficiency compared to the conventional linear state transition matrix with Euler or Runge-Kutta integration. © 2010 The Japan Society for Aeronautical and Space Sciences.