  In this paper we apply the Regge-Wheeler formalism to study the propagation
of axial and polar gravitational waves in matter-filled Bianchi I universe.
Assuming that the expansion scalar $ \Theta $, of the background space-time, is
proportional to the shear scalar $ \sigma $, we solved the background field
equations in the presence of matter (found to behave like a stiff fluid). We
then derive the linearised perturbation equations for both the axial and polar
modes. The analytical solutions in vacuum spacetime could be determined in an
earlier paper \cite{GD1} in a relatively straightforward manner. However, here
we find that in the presence of matter, they require more assumptions for their
solution, and bear more involved forms. As compared to the axial modes, the
polar perturbation equations contain far more complicated couplings among the
perturbing terms. Thus we have to apply suitable assumptions to derive the
analytical solutions for some of the cases of polar perturbations. In both the
axial and polar cases, the radial and temporal solutions for the perturbations
separate out as products. We find that the axial waves are damped owing to the
background anisotropy, and can deform only the azimuthal velocity of the fluid.
In contrast, the polar waves must trigger perturbations in the energy density,
the pressure as well as in the non-azimuthal components of the fluid velocity.
Similar behaviour is exhibited by axial and polar gravitational waves
propagating in the Kantowski-Sachs universe \cite{GD2}. Our work is in contrast
to the work done in \cite{SYK}, where the authors analysed anisotropic
universes modelled by Kasner spacetime and Rindler wedges using the method of
gauge-invariant perturbations in the RW gauge.
