In recent years, Com–Poisson has emerged as one of the most popular discrete models in the analysis of count data owing to its flexibility in handling different types of dispersion. However, in a stationary longitudinal Com–Poisson count data set-up where the covariates are time independent, estimation of regression and dispersion parameters based on a generalized quasi-likelihood (GQL) approach involves some major computational difficulties particularly in the inversion of the joint covariance matrix. On the other hand, in practical real-life longitudinal studies, time-dependent covariates leading to non-stationary responses are more frequently encountered. This implies that further computational problems will now arise when estimating parameters under non-stationary set-ups. This paper overcomes this problem by approximating the inverse of the ill-conditioned covariance matrix in the GQL approach through a multidimensional conjugate gradient method. The performance of this novel version of the GQL approach is then assessed on simulations of AR(1) stationary and AR(1) non-stationary longitudinal Com–Poisson counts and on real-life epileptic seizure counts. However, there is not yet an algorithm to generate non-stationary longitudinal Com–Poisson counts nor a GQL algorithm to estimate the parameters under non-stationary set-ups. Thus, the paper also provides a framework to generate non-stationary AR(1) Com–Poisson counts along with the construction of a GQL equation under non-stationary set-ups.