We build upon recent advances on the distributional aspect of Stein's method to propose a novel and flexible technique for computing Stein operators for random variables that can be written as products of independent random variables. We show that our results are valid for a wide class of distributions including normal, beta, variance-gamma, generalized gamma and many more. Our operators are kth degree differential operators with polynomial coefficients; they are straightforward to obtain even when the target density bears no explicit handle. As an application, we derive a new formula for the density of the product of k independent symmetric variance-gamma distributed random variables. -Tilman, Allée de la découverte 12, B-8000 Liège, Belgium yswan@uliege.be. variable of interest W and taking the supremum over all h in some class of functions H leads to the estimatewhere the final supremum is taken over all f h that solve (1). The third and final step of the method involves developing appropriate strategies for bounding the expectation on the right hand side of (2). This is of interest because many important probability metrics (such as the Kolmogorov and Wasserstein metrics) are of the form d H (L(W ), L(X)). Moreover, in many settings bounding the expectation IE[Af h (W )] is relatively tractable, and as a result Stein's method has found application in disciplines as diverse as random graph theory [5], number theory [22], statistical mechanics [13] and quantum mechanics [29]. We refer to the survey paper [37] as well as to the monographs [30,9] for a deeper look into some of the fruits of Charles Stein's seminal insights, particularly in the case where the target is the normal distribution. The linchpin of the method is the operator A whose properties are crucial to the success of the whole enterprise. In the sequel, we concentrate exclusively on differential Stein operators (some operators in the literature are integral or even fractional, see e.g. [43,3]) and adopt the following lax definition:Definition 1.1. A linear differential operator A acting on a class F of functions is a Stein operator for X if (i) Af ∈ L 1 (X) and (ii) IE [Af (X)] = 0 for all f ∈ F.There are infinitely many Stein operators for any given target distribution. For instance, if the distribution is known (even if only up to a normalizing constant) then the "canonical" theory from [26] applies, leading to entire families of operators. This approach provides natural first order polynomial operators e.g. for target distributions which belong to the Pearson family [38] or which satisfy a diffusive assumption [11,25]. In some cases, one may rather apply a duality argument. For instance the p.d.f. γ(x) = (2π) −1/2 e −x 2 /2 of the standard normal distribution satisfies the first order ODE γ ′ (x) + xγ(x) = 0 leading, by integration by parts, to the already mentioned operator Af (x) = f ′ (x) − xf (x). This is particularly useful for densities defined implicitely via ODEs. Such are by no means the only methods for deriving differential Stein ope...