  We establish the asymptotic expansion in $\beta$ matrix models with a
confining, off-critical potential, in the regime where the support of the
equilibrium measure is a union of segments. We first address the case where the
filling fractions of these segments are fixed, and show the existence of a
$1/N$ expansion. We then study the asymptotics of the sum over the filling
fractions, to obtain the full asymptotic expansion for the initial problem in
the multi-cut regime. In particular, we identify the fluctuations of the linear
statistics and show that they are approximated in law by the sum of a Gaussian
random variable and an independent Gaussian discrete random variable with
oscillating center. Fluctuations of filling fractions are also described by an
oscillating discrete Gaussian random variable. We apply our results to study
the all-order small dispersion asymptotics of solutions of the Toda chain
associated with the one Hermitian matrix model ($\beta = 2$) as well as
orthogonal ($\beta = 1$) and skew-orthogonal ($\beta = 4$) polynomials outside
the bulk.
