How to average coherences estimated using Welch's method in MATLAB
If you look at the documentation of mscohere which contains the definition of magnitude squared coherence, you will see that it is the absolute squared of the cross spectral density, divided by the product of the two separate spectral densities. Since power densities enter numerator and denominator of the expression, the MSC estimate from a full time series is not simply the average of the estimates from the parts.

To get a good "average" MSC estimate, do the following: estimate in each segment separately cross and single spectral densities, average these three over the segments, and then combine them in the end into an average MSC. Continuing your code:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);
sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);
sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;
sxya = (sxy1 + sxy2) / 2;
ca = abs(sxya) .^ 2 ./ sxa ./ sya;


Now max(abs(c-ca)) is almost zero.