While DNA sequence evolution has been well studied, the expression of genes is also subject to evolution. Yet the evolution of gene expression is currently not well understood. In recent years, new tissue/organ specific gene expression datasets spanning several organisms across the tree of life, have become available providing the opportunity to study gene expression evolution in more detail. However, while a theoretical model to study evolution of continuous traits exist, in practice computational methods often cannot distinguish, with confidence, between alternative evolutionary scenarios. This lack of power has been attributed to the modest number of species with available expression data.To solve this challenge, we introduce EvoGeneX, a computationally efficient method to uncover the mode of gene expression evolution based on the Ornstein-Uhlenbeck process. Importantly, EvoGeneX in addition to modelling expression variations between species, models within species variation. To estimate the within species variation, EvoGeneX formally incorporates the data from biological replicates as a part of the mathematical model. We show that by modelling the within species diversity EvoGeneX significantly outperforms the currently available computational method. In addition, to facilitate comparative analysis of gene expression evolution, we introduce a new approach to measure the dynamics of evolutionary divergence of a group of genes.We used EvoGeneX to analyse the evolution of expression across different organs, species and sexes of the Drosophila genus. Our analysis revealed differences in the evolutionary dynamics of male and female gonads, and uncovered examples of adaptive evolution of genes expressed in the head and in the thorax.