  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 [Phys. Rev. Lett. 126, 021602]. 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 [Phys. Rev.
Lett. 126, 021602]. For 5-point conformal blocks, there exist five such
operators which are worked out smoothly in the dimension $d$.
