A numerical method is developed for solving the two-dimensional Helmholtz equation in a two-layer region bounded by a flat top, a flat bottom and a curved interface. A better local orthogonal transformation is used to flatten the curved interface of the waveguide. The one-way re-formulation based on the Dirichlet-to-Neumann (DtN) map is then used to reduce the boundary value problem to an initial value problem. Numerical implemetation of the resulting operator Riccati equation uses a large range step method for discretizing the range variable and a truncated local eigenfunction expansion for approximating the operators. This method is particularly useful for solving long range wave propagation problems in slowly varying waveguides.