The parabolic equation technique provides an excellent model to describe the wave phenomena when there exists a predominant direction of propagation. The model handles the square root wave number operator in paraxial direction. Realization of the pseudo-differential square root operator is the essential part of the parabolic equation method for its numerical accuracy. The wide-angled approximation of the operator is made based on the Pade series expansion, where the branch line rotation scheme can be combined with the original Pade approximation to stabilize its computational performance for complex modes. The Galerkin integration has been employed to discretize the depth-dependent operator. The benchmark tests involving the half-infinite space, the range independent and dependent environment will validate the implemented numerical model.