Fractional Brownian motion is a Gaussian stochastic process with stationary, long-time correlated increments and is frequently used to model anomalous diffusion processes. We study numerically fractional Brownian motion confined to a finite interval with reflecting boundary conditions. The probability density function of this reflected fractional Brownian motion at long times converges to a stationary distribution showing distinct deviations from the fully flat distribution of amplitude $1/L$ in an interval of length $L$ found for reflected normal Brownian motion. While for superdiffusion, corresponding to a mean squared displacement $langle X^2(t)ranglesimeq t^{alpha}$ with $1<alpha<2$, the probability density function is lowered in the centre of the interval and rises towards the boundaries, for subdiffusion ($0<alpha<1$) this behaviour is reversed and the particle density is depleted close to the boundaries. The mean squared displacement in these cases at long times converges to a stationary value, which is, remarkably, monotonically increasing with the anomalous diffusion exponent $alpha$. Our a priori surprising results may have interesting consequences for the application of fractional Brownian motion for processes such as molecule or tracer diffusion in the confined of living biological cells or organelles, or other viscoelastic environments such as dense liquids in microfluidic chambers.