Transmission of natural gas through high pressure pipelines has been modeled by numerically solving the governing equations for one-dimensional compressible flow using implicit finite difference methods. In the first case the backward Euler method is considered using both standard first-order upwind and second-order centered differences for the spatial derivatives. The first-order upwind approximation, which is a one-sided approximation, is found to be unstable for CFL numbers less than $1$, while the centered difference approximation is stable for any CFL number. In the second case a cell centered method is considered where flow values are calculated at the midpoint between grid points. This method is also stable for any CFL number. However, for a discontinuous change in inlet temperature, the method is observed to introduce unphysical oscillations in the temperature profile along the pipeline. A solution strategy where the hydraulic and thermal models are solved separately using different discretization techniques is suggested. Such a solution strategy does not introduce unphysical oscillations for discontinuous changes in inlet boundary conditions and is found to be stable for any CFL number. The one-dimensional flow model is validated using operational data from a high pressure natural gas pipeline.