We present a numerical study of the propagation of VLF whistler waves in the magnetospheric plasma. In this study the plasma is considered to be homogeneous in the direction along the ambient magnetic field and strongly inhomogeneous across it. The goal of this investigation is to understand whistler propagation in magnetic-field-aligned channels (also called ducts) with either enhanced or depleted plasma density. In particular, the paper is focused on situations where the transverse scale size of the duct is comparable to or smaller than the perpendicular wavelength of the whistler. In this case, classical analysis of the whistler dynamics based on the geometrical optic approximation becomes invalid, and numerical solutions of the full wave equations should be performed. Our simulations extend the earlier analysis based on the ray-tracing technique and analytical studies of the very low frequency wave equations. We show that high-density ducts are inherently leaky and this leakage depends on the perpendicular wavelength of the wave inside the duct. We also show that whistler trapping occurs not only at density maxima and minima but also at critical points along a density gradient. This effect can explain whistler guiding along strong transverse plasma density gradients at the plasmapause.