The heliopause is the interface between the solar wind plasma and the very local interstellar medium (VLISM) and manifests itself as a tangential discontinuity across which the flow velocity and the plasma density jump (except at the nose). Hydrodynamic instabilities of either the Rayleigh‐Taylor type or the Kelvin‐Helmholtz type will likely develop at the heliopause. To our knowledge, previous analytical studies of these instabilities were confined to linear perturbation analyses, and most existing numerical simulations did not obtain the Kelvin‐Helmholtz type instability of the heliopause, probably due to large numerical dissipation. In this paper we use the piecewise parabolic method (PPM) in our hydrodynamic simulation to study the stability of the heliopause. The PPM can capture shocks and discontinuities within 1–2 grid points with negligible numerical dissipation. For simplicity, magnetic fields, interstellar neutrals, cosmic rays, etc., are neglected in our model. We thus focus our attention on the general pattern of the Kelvin‐Helmholtz instability at the heliopause. In both the “one‐shock” and “two‐shock” models, the Kelvin‐Helmholtz instability occurs at the heliopause and leads to nonlinear oscillations of the heliopause and the termination shock with a timescale of the order of 102 years. The excursion of the heliopause at the nose as a result of these oscillations is of the order of tens of astronomical units, with much smaller excursions for the termination shock. Growth rates from the simulations are in reasonable agreement with theoretical estimates. The possible stabilizing influence of the magnetic field, neglected in the present model, is discussed.