The first essential ingredient to build up Stein's method for a continuous target distribution is to identify a so-called Stein operator, namely a linear differential operator with polynomial coefficients. In this paper, we introduce the notion of algebraic Stein operators (see Definition 3.2), and provide a novel algebraic method to find all the algebraic Stein operators up to a given order and polynomial degree for a target random variable of the form Y = h(X), where X = (X 1 , . . . , X d ) has i.i.d. standard Gaussian components and h ∈ K[X] is a polynomial with coefficients in the ring K. Our approach links the existence of an algebraic Stein operator with null controllability of a certain linear discrete system. A MATLAB code checks the null controllability up to a given finite time T (the order of the differential operator), and provides all null control sequences (polynomial coefficients of the differential operator) up to a given maximum degree m. This is the first paper that connects Stein's method with computational algebra to find Stein operators for highly complex probability distributions, such as H 20 (X 1 ), where H p is the p-th Hermite polynomial. A number of examples of Stein operators for H p (X 1 ), p = 3, 4, 5, 6, 7, 8, 10, 12, are gathered in the extended Appendix of this arXiv version. We also introduce a widely applicable approach to proving that Stein operators characterise the target distribution, and use it to prove, amongst other examples, that the Stein operators for H p (X 1 ), p = 3, . . . , 8, with minimum possible maximal polynomial degree m characterise their target distribution.