The sound propagation problems in range-dependent waveguides are a common topic in underwater acoustics. The range-dependent factors, involving volumetric and bathymetric variations, significantly influence the propagation of sound energy and information. In this paper, a coupled-mode method based on the multimodal admittance method is presented for analyzing the sound propagation and scattering problems in range-dependent waveguides. The sound field is expanded in terms of a local basis with range-dependent modal amplitudes. The local basis corresponds to the transverse modes in a waveguide with constant physical parameters and constant cross section equal to the local cross section in the range-dependent waveguide. This local basis takes the advantage that it is easier to compute than the usual local modes which are the transverse modes in a waveguide with local physical parameters and constant cross-section equal to the local cross-section, especially for waveguides with complex environments. Projection of the Helmholtz equation that governs the sound pressure onto the local basis gives the second-order coupled mode equations for the modal amplitudes of the sound pressure. The correct boundary conditions are used in the derivation, giving rising to boundary matrices, in order to guarantee the conservation of energy among modes. The second-order coupled mode equations include coupled matrices and boundary matrices, which directly describe the effect of mode coupling due to contribution from volumetric variation (range-dependent physical parameters) and bathymetric variation (range-dependent boundaries). By introducing the admittance matrix, the second-order coupled mode equations are reduced to two sets of first-order evolution equations. The Magnus integration method is used to solve the first-order evolution equations. These first-order evolution equations allow us to obtain the numerical stable solutions and avoid the numerical divergence due to the exponential growth of evanescent modes. The numerical examples are presented for the waveguides with range-dependent physical parameters or range-dependent boundaries. The agreement between the results computed with the coupled mode method and COMSOL verifies the accuracy of the coupled mode method. Although the analysis and numerical implementation in this paper are based on two-dimensional waveguides in Cartesian coordinate system, it can be generally extended to study more complex waveguides.