We present a novel numerical model that simulates the anisotropic Bragg diffraction in optically anisotropic (uniaxial) acousto-optical devices. We use a non-paraxial vectorial beam propagation method adapted to optically inhomogeneous medium and arbitrary optical field distribution. The principal idea of our solving method is that since the amplitude of the spatial variation of the refractive index caused by the acoustic wave is relatively small, we can consider it as a perturbation and iterate to the exact solution of the wave equation describing the propagation of the optical field distribution. To describe anisotropic diffraction with polarization rotation, we use a fully vectorial beam propagation method (BPM) where the accuracy depends on the relative step size. The results converge rapidly to fulfill energy conservation (up to 10-3 with less than an hour of computational time). We show that the calculated angles, the space dependent intensity and polarization variations of the diffracted beams agree accurately with those predicted by theory and experiment, under Bragg diffraction condition, in various acousto-optic configurations.