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 5dimensional 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+δ , with δ 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.