A novel vectorial finite-element beam propagation method for anisotropic medium is presented. The approach has been derived in order to be suitable for the numerical discretization through the finite elements and able to take into account both polarization dependence and field component coupling has been derived. The vector wave equation is solved in terms of the transverse magnetic field components, such that the spurious modes are eliminated by an implicit inclusion of divergence condition. The precision of the proposed approach is demonstrated through numerical examples.