The authors present an automated approach to constructing outer approximations for systems in a class of integro-difference operators with smooth nonlinearities. Chebyshev interpolants and Galerkin projections form the basis for the construction, while analysis and interval arithmetic are used to incorporate explicit error bounds. The results represent a significant advance compared to the approach given by S. Day et al. [SIAM J. Appl. Dyn. Syst. 3, No. 2, 117–160 (2004; Zbl 1059.37068)], extending the nonlinearities that may be studied from low-degree polynomials to smooth functions and the studied portion of phase space from a simulated attracting region to the global maximal invariant set.