Co-existence of the physical and numerical boundary conditions makes implicit boundary treatment a particularly difficult problem in modern CFD simulations. Previous studies adopted space–time extrapolation or specially designed partial differential equations on the boundaries that are different from those of interior points. They are often formulated in terms of primitive variables, and are very challenging for complicated boundary types to be converted to, and implemented in, the conservative variables that are preferred in numerical simulations. More importantly, different boundary equations or different extrapolation techniques may compromise the stability, accuracy, or convergence rate of the A-stable schemes that are developed for interior points. A new methodology for implicit boundary treatment is proposed in this study. By introducing a simple correction matrix T, a set of generalized equations ∂Q/∂t=(I+T)R are developed in terms of conservative variables. It is applicable for both the interior domain (T=0) and the boundaries (T≠0). It is in a partial differential equation form, satisfies the boundary conditions accurately, independent of the time and spatial discretizations. Any one-sided schemes can be used on the boundaries but still maintain the upwind property. Implicit solution techniques are made significantly easy to implement using, for example, the data-parallel lower–upper relation method and the Newton method (combined with the GMRES method for subsidiary iterations). Numerical experiments show that the proposed methodology produces stable simulations for very large CFL numbers and preserve the imposed boundary values accurately.