Like Damian has said, the BP method is suitable for propagation of length up to several thousand times the wavelength. (So for optical wavelengths, millimetre lengths of propagation). For longer propagation, you’ll need to use other methods. There is probably no commercial software for that simulation, but you can probably make a suitable model yourself. For this problem, I recommend non-orthogonal coupled mode theory. When solving by that method, you will need to know details of the optical properties of the waveguide, such as modes, modal index, modal field pattern. You’ll need to calculate overlap integrals of the optical modal fields with the coupler geometry. Optiwave software can help by finding the modes and also with overlap integral calculation. We don’t have the complete solution for you, but our software can help with the hardest part: finding the modes.
.
.