Bayesian inference is conceptually elegant, but calculating posterior expectations can entail a heavy computational cost. Monte Carlo methods are reliable and supported by strong asymptotic guarantees, but do not leverage smoothness of the integrand. Solving Stein equations has emerged as a possible alternative, providing a framework for numerical approximation of posterior expectations in which smoothness can be exploited. However, existing numerical methods for Stein equations are associated with high computational cost due to the need to solve large linear systems. This paper considers the combination of iterative linear solvers and preconditioning strategies to obtain fast approximate solutions of Stein equations.