Abstract. A computational method based on Chebyshev series to rigorously compute solutions of initial and boundary value problems of analytic nonlinear vector fields is proposed. The idea is to recast solutions as fixed points of an operator defined on a Banach space of rapidly decaying Chebyshev coefficients and to use the so-called radii polynomials to show the existence of a unique fixed point nearby an approximate solution. As applications, solutions of initial value problems in the Lorenz equations and symmetric connecting orbits in the Gray-Scott equation are rigorously computed. The symmetric connecting orbits are obtained by solving a boundary value problem with one of the boundary values in the stable manifold.Key words. Rigorous numerics, chebyshev series, nonlinear ODEs, boundary value problems, initial value problems, contraction mapping theorem AMS subject classifications. 34B99, 34A34, 34A12, 42A10, 65L05, 65L101. Introduction. In this paper, we propose a rigorous numerical method based on the Chebyshev polynomials to compute solutions of nonlinear differential equations. More explicitly, the field of rigorous numerics develops algorithms that provide approximate solutions to a problem together with precise bounds within which exact solutions are guaranteed to exist in the mathematically rigorous sense. In this context, the main idea of our proposed approach is to expand the solution of a given differential equation using its Chebyshev series, plug the expansion in the equation, obtain an equivalent infinite dimensional problem of the form f (x) = 0 to solve in a Banach space of rapidly decaying Chebyshev coefficients and to get the existence, via a fixed point argument, of a genuine solution of f (x) = 0 nearby a numerical approximation of a finite dimensional projection of f . The fixed point argument is solved by using the radii polynomials (e.g. see [1]), which provide an efficient way of constructing a set on which the contraction mapping theorem is applicable.Before proceeding further, it is worth mentioning that a similar approach based on Fourier series is widely used in the field of rigorous numerics to compute solutions of differential equations with periodic profiles. For instance, time periodic solutions of ODEs [2,3], stationary solutions of PDEs with periodic or Neumann boundary conditions [4,5,6,7], time periodic solutions of delay differential equations [8,9] and invariant sets of infinite dimensional maps [10] have been successfully computed using Fourier series and rigorous numerics. However, to the best of our knowledge, this is the first time that a method based on Chebyshev series is presented to rigorously compute solutions of nonlinear differential equations. Since a large class of solutions of differential equations are non periodic (e.g. solutions of initial value problems (IVPs) and boundary value problems (BVPs) with non periodic boundary values), we believe that our proposed approach is a valuable contribution to the field of rigorous numerics. Also, since Chebyshev seri...