We have determined the thermal expansion of the first three interlayer spacings of Be(0001) by comparing measured temperature-dependent surface core level shifts and surface state positions to values calculated with Density Functional Theory for different surface geometries. Our results reveal that, in the temperature range from 300 to 700 K, the 1st-to-2nd, 2nd-to-3rd and 3rd-to-4th interlayer coefficient of thermal expansion are 88 +/- 15, -10 +/- 15 and -6 +/- 20 x 10(-6) K(-1), respectively. We have now set Be(0001) as a challenge for further theoretical improvements, since a complete theoretical understanding of surface thermal expansion does not exist.