An implicit finite difference method is described for solving the equations of motion of a two-layer stably stratified fluid in one dimension. The computational scheme is applicable to smoothly varying subcritical flows and is formulated in terms of layer transport and thickness in order to optimize computer efficiency. Stability analyses of the method show that waves forming at the grid scale are undamped and two approaches for controlling this wave growth are discussed. Comparisons with experimental data show the method to be accurate provided that the natural long wave motions are well resolved and the grid scale waves are filtered from the solutions. Numerical tests indicate a minimum of 20 points per wave length or period can be considered as a lower limit for resolution.