We present a gradient-based topology optimization method of optical waveguide devices using the beam propagation method (BPM). Efficient topology optimization of an optical waveguide is available by using our approach because it utilizes the adjoint variable method for sensitivity analysis, and the BPM for numerical propagation analysis. In this paper, a way of evaluating sensitivities for a design of longitudinally varying waveguides is newly described in the case that a density method and finite difference BPM are employed. In addition, the validity and availability of our proposed method are verified by optimizing an S-bend and a Y-branch waveguides.