Many evolutionary processes can lead to a change in the correlation between continuous characters over time or on different branches of a phylogenetic tree. Shifts in genetic or functional constraint, in the selective regime, or in some combination thereof can influence both the evolution of continuous traits and their relation to each other. These changes can often be mapped on a phylogenetic tree to examine their influence on multivariate phenotypic diversification. We propose a new likelihood method to fit multiple evolutionary rate matrices (also called evolutionary variance-covariance matrices) to species data for two or more continuous characters and a phylogeny. The evolutionary rate matrix is a matrix containing the evolutionary rates for individual characters on its diagonal, and the covariances between characters (of which the evolutionary correlations are a function) elsewhere. To illustrate our approach, we apply the method to an empirical dataset consisting of two features of feeding morphology sampled from 28 centrarchid fish species, as well as to data generated via phylogenetic numerical simulations. We find that the method has appropriate type I error, power, and parameter estimation. The approach presented herein is the first to allow for the explicit testing of how and when the evolutionary covariances between characters have changed in the history of a group.