Seismic modeling plays an important role in geophysics and seismology for estimating the response of seismic sources in a given medium. In this work, we present a MATLAB-based package, FDwave3D, for synthetic wavefield and seismogram modeling in 3D anisotropic media. The seismic simulation is carried out using the finite-difference method over the staggered grid, and it is applicable to both active and passive surveys. The code package allows the incorporation of arbitrary source mechanisms and offers spatial derivative operators of accuracy up to tenth-order along with different types of boundary conditions. First, the methodological aspects of finite-difference method are briefly introduced. Then, the code has been tested and verified against the analytical solutions obtained for the homogeneous model. Further, the numerical examples of layered and overthrust models are presented to demonstrate its reliability.