An idea which has been around in general relativity for more than forty years is that in the approach to a big bang singularity solutions of the Einstein equations can be approximated by the Kasner map, which describes a succession of Kasner epochs. This is already a highly non-trivial statement in the spatially homogeneous case. There the Einstein equations reduce to ordinary differential equations and it becomes a statement that the solutions of the Einstein equations can be approximated by heteroclinic chains of the corresponding dynamical system. For a long time progress on proving a statement of this kind rigorously was very slow but recently there has been new progress in this area, particularly in the case of the vacuum Einstein equations. In this paper we generalize some of these results to the Einstein-Maxwell equations. It turns out that this requires new techniques since certain eigenvalues are in a less favourable configuration in the case with a magnetic field. The difficulties which arise in that case are overcome by using the fact that the dynamical system of interest is of geometrical origin and thus has useful invariant manifolds.