Existing algorithms for bidirectional optical beam propagation proposed to simulate reflective integrated photonic devices do not propagate evanescent fields correctly. Thus inaccuracy and instability problems can arise when fields have significant evanescent character. We propose complex representations of the propagation operator by choosing either a complex reference wave number or a complex representation of Pade approximation to address this issue. Therefore correct evolution of both propagating waves and evanescent waves can be simultaneously realized, significantly reducing the inaccuracy and instability problems. Both test problems and practical problems are presented for demonstration.