We have developed a population density framework that captures correlations between any pair of neurons in the population. We model each population of integrate-and-fire neurons as receiving input in the form of correlated Poisson processes. The evolution equation for the probability density of any pair of neurons within the population is a multivariate integro-differential equation which we solve numerically. We demonstrate the numerical method and compare the numerical solutions with Monte-Carlo simulations. Traditional population density approaches assume all neurons within a population are independent. However, correlations that are missed by these approaches can significantly alter network dynamics. Hence, the correlated population density method developed here could provide a framework to analyze how correlations propagate through networks and could be a computationally efficient method to accurately simulate large scale networks.