We introduce a simple iterative technique for bounding derivatives of solutions of Stein equations $Lf=h-\mathbb{E}h(Z)$, where $L$ is a linear differential operator and $Z$ is the limit random variable. Given bounds on just the solutions or certain lower order derivatives of the solution, the technique allows one to deduce bounds for derivatives of any order, in terms of supremum norms of derivatives of the test function $h$. This approach can be readily applied to many Stein equations from the literature. We consider a number of applications; in particular, we derive bounds for derivatives of any order of the solution of the general variance-gamma Stein equation.