The present work analyzes stationary distributions of active Brownian particles in a harmonic trap. Generally, obtaining stationary distributions for this system is non-trivial, and to date, no exact expressions are available. In this work, we develop and explore a method based on a transformation of the Fokker–Planck equation into a recurrence relation for generating moments of a distribution. The method, therefore, offers an analytically tractable approach, an alternative to numerical simulations, in a situation where more direct analytical approaches fail. Although the current work focuses on the active Brownian particle model, the method is general and valid for any type of active dynamics and any system dimension.