Abstract The stochastic solution to diffusion equations with polynomial coefficients is called a Pearson diffusion. If the time derivative is replaced by a distributed order fractional derivative, the stochastic solution is called a distributed order fractional Pearson diffusion. This paper develops a formula for the covariance function of distributed order fractional Pearson diffusion in the steady state, in terms of generalized Mittag-Leffler functions. The correlation function decays like a power law. That formula shows that distributed order fractional Pearson diffusions exhibits long range dependence.