We present a method based on the Grüneisen formalism to calculate the thermal expansion coefficient (TEC) tensor that is applicable to any crystal system, where the number of phonon calculations associated with different deformations scales linearly with the number of lattice parameters. Compared to simple high-symmetry systems such as cubic or hexagonal systems, a proper consideration of low-symmetry systems such as monoclinic or triclinic crystals demands a clear distinction between the TEC tensor and the lattice-parameter TECs along the crystallographic direction. The latter is more complicated and it involves integrating the equations of motion for the primitive lattice vectors, with input from the TEC tensor. A first-principles study of the TEC is carried out for the first time on a monoclinic crystal, where we unveil high TEC anisotropies in a recently reported monoclinic phase of niobium trisulfide (NbS3-IV) crystal with a relatively large primitive cell (32 atoms per cell) using density-functional theory. We find the occurrence of a negative TEC tensor component is largely due to the mechanical property rather than the anharmonic effect, contrary to common belief. Our theoretical treatment of the monoclinic system with a single off-diagonal tensor element could be routinely generalized to any crystal system, including the lowest-symmetry triclinic system with three off-diagonal tensor elements.