Multiple and co-evolving networks are common in many real settings such as social networks, communication networks and other information networks. Most of the work in the field of network evolution has focused on a single evolving network or specific network pairs, lacking generality in the analysis of multiple networks and ignoring the co-evolutionary dynamics between networks. In practice, a significant amount of information is encoded in the evolution of multiple networks with respect to one another. In this paper, we show how to use a shared temporal matrix factorization framework to model co-evolution across multiple networks, and we refer to this framework as CoEvOL. Specifically, the proposed framework decomposes the adjacency matrix of each co-evolving network into a product of network-independent shared factor and a set of network-specific temporal factors, and impose a non-negativity constraint on the factors for greater interpretability. Our approach has the potential to predict multiple changes in co-evolving networks over time, because of its ability to explicitly represent co-evolving networks as a function of time. The CoEvoL framework also has the advantage of generality in addressing various temporal tasks across multiple networks. We show the benefits of this approach in predicting co-evolution across multiple networks on the tasks including cross-network link prediction, lag correlation detection and community detection. Compared to baseline methods, CoEvoL obtains lower root mean-squared error in cross-network link prediction and higher cluster purity in community detection, which demonstrates that the Co-EvoL framework can capture the dynamics across multiple networks.