  We present a dynamical analysis in terms of new expansion-normalized
variables for homogeneous and anisotropic Bianchi-I spacetimes in $f(R)$
gravity in the presence of anisotropic matter. With a suitable choice of the
evolution parameter, the Einstein's equations are reduced to an autonomous
5-dimensional system of ordinary differential equations for the new variables.
Further restrictions lead to considerable simplifications. For instance, we
show that for a large class of functions $f(R)$, which includes several cases
commonly considered in the literature, all the fixed points are polynomial
roots, and hence they can be determined with good accuracy and classified for
stability. Moreover, typically for these cases, any fixed point corresponding
to isotropic solutions in the presence of anisotropic matter will be unstable.
The assumption of a perfect fluid as source and or the vacuum cases imply some
dimensional reductions and even more simplifications. In particular, we find
that the vacuum solutions of $f(R) = R^{1+\delta}$, with $\delta$ a constant,
are governed by an effective bi-dimensional phase space which can be
analytically constructed, leading to an exactly soluble dynamics. Finally, we
demonstrate that several results already reported in the literature can be
re-obtained in a more direct and easy way by exploring our dynamical
formulation.
