A 3D full-vectorial Beam Propagation Method is successfully applied to compute both the propagation constants and the modal profiles in high-contrast silica-air index-guiding Photonic Crystal Fibers. The approach is intrinsically suited to investigate longitudinally varying structures or propagation and polarization effects, which are of practical interest for advanced optical applications. As an example we model a dual-core coupler, showing that efficient polarization preserving coupling can be expected.