The construction of conformal blocks for the analysis of multipoint correlation functions with N > 4 local field insertions is an important open problem in higher dimensional conformal field theory. This is the first in a series of papers in which we address this challenge, following and extending our short announcement in [1]. According to Dolan and Osborn, conformal blocks can be determined from the set of differential eigenvalue equations that they satisfy. We construct a complete set of commuting differential operators that characterize multipoint conformal blocks for any number N of points in any dimension and for any choice of OPE channel through the relation with Gaudin integrable models we uncovered in [1]. For 5-point conformal blocks, there exist five such operators which are worked out smoothly in the dimension d.