We present results of numerical simulations that examine the dynamical stability of known planetary systems, a star with two or more planets. First we vary the initial conditions of each system on the basis of observational data. We then determine regions of phase space that produce stable planetary configurations. For each system we perform 1000∼106 yr integrations. We examine υ And, HD 83443, GJ 876, HD 82943, 47 UMa, HD 168443, and the solar system. We find that the resonant systems, two planets in a first-order mean motion resonance (HD 82943 and GJ 876) have very narrow zones of stability. The interacting systems, not in first-order resonance, but able to perturb each other (υ And, 47 UMa, and the solar system), have broad stable regions. The separated systems, two planets beyond 10:1 resonance (we examine only HD 83443 and HD 168443) are fully stable. We find that the best fits to the interacting and resonant systems place them very close to unstable regions. The boundary in phase space between stability and instability depends strongly on the eccentricities and (if applicable) the proximity of the system to perfect resonance. Furthermore, we also find that the longitudes of periastron circulate in chaotic systems but librate in regular systems. In addition to 106 yr integrations, we also examined stability on ∼108 yr timescales. For each system we ran ∼10 long-term simulations, and find that the Keplerian fits to these systems all contain configurations that are regular on this timescale.