Fast and accurate simulation of responses of electromagnetic logging-while-drilling tools in 2D and 3D formations is important for further development of the proactive geosteering service, which is currently based on the 1D layered-model inversion. In search of an optimal method of forward simulation, we have developed boundary integral equations in a 2D formation model with planar boundaries. A peculiarity of the method is that the number of grid nodes required for reaching a certain accuracy, and hence the computation time increases as the transmitter approaches the boundary. To retain efficiency in the case of transmitter-boundary proximity, we use the following modification of the method. The main term of the calculated electromagnetic field, namely, the solution for the two-layered model corresponding to the closest boundary, is calculated by the well-known explicit formulas, whereas the anomalous field is found from the system of integral equations. The modified method is implemented and tested on complex realistic 2D formation models as well as basic models to be used in 2D inversion. Numerical simulations of real propagation and azimuthal resistivity tools show that software based on the modified method can be successfully used for synthetic log simulation in complex formations. However, the proposed method has two shortcomings: poor efficiency when the transmitter is close to the boundary endpoints and limitation to isotropic formations. Both are topics for further research. Currently, the singularity near the endpoints is treated by simply refining the grid, which deteriorates efficiency. Generalization to the anisotropy case would include modification of the basic method and derivation of the closed-form solution for the layered model without the assumption of transverse isotropy.