Proving the bound: The sample variance is not affected by location shifts in the data, so we will use the values $0 \leqslant y_i \leqslant b-a$ and compute the sample variance from these values. As a preliminary observation we first note that the bound $y_i \leqslant b-a$ implies that we have $\sum y_i^2
\leqslant \sum (b-a) y_i = (b-a) n \bar{y}_n$, which gives the upper bound:
$$\begin{align}
s_n^2
&= \frac{1}{n-1} \sum_{i=1}^n (y_i - \bar{y}_n)^2 \\[6pt]
&= \frac{1}{n-1} \Bigg[ \sum_{i=1}^n y_i^2 - n \bar{y}_n^2 \Bigg] \\[6pt]
&\leqslant \frac{1}{n-1} \Bigg[ (b-a) n \bar{y}_n - n \bar{y}_n^2 \Bigg] \\[6pt]
&= \frac{n}{n-1} \cdot \bar{y}_n (b-a - \bar{y}_n). \\[6pt]
\end{align}$$
Now, using the fact that $0 \leqslant \bar{y}_n \leqslant b-a$ this upper bound is maximised when $\bar{y}_n = \tfrac{b-a}{2}$, which gives:
$$\begin{align}
s_n^2
&\leqslant \frac{n}{n-1} \cdot \bar{y}_n (b-a - \bar{y}_n) \quad \ \ \ \\[6pt]
&\leqslant \frac{n}{n-1} \cdot \Big( \frac{b-a}{2} \Big)^2. \\[6pt]
\end{align}$$
This gives an upper bound for the sample variance. Splitting the points between the two bounds $a$ and $b$ (for an even number of points) achieves this upper bound, so it is the maximum possible sample variance for bounded data.
Maximum possible sample variance with an odd number of data points: It is also possible to find the upper bound when we have an odd number of data points. One way to do this is to use the iterative updating formula for the sample variance (see e.g., O'Neill 2014). If we consider the maximising case for an even number $n-1$ of values and then take the final value to be $x_n = x$, then we have:
$$\begin{align}
\quad \quad \quad \quad \quad \quad s_n^2
&= \frac{1}{n-1} \Bigg[ (n-2) s_{n-1} + \frac{n-1}{n} \cdot (x-\bar{x}_{n-1})^2 \Bigg] \\[6pt]
&\leqslant \frac{1}{n-1} \Bigg[ (n-1) \cdot \Big( \frac{b-a}{2} \Big)^2 + \frac{n-1}{n} \cdot \Big( x-\frac{b-a}{2} \Big)^2 \Bigg] \\[6pt]
&= \Big( \frac{b-a}{2} \Big)^2 + \frac{1}{n} \cdot \Big( x-\frac{b-a}{2} \Big)^2. \\[6pt]
\end{align}$$
This quantity is maximised by taking either $x=a$ or $x=b$ (i.e., the last point is also on the boundary), which gives:
$$\begin{align}
s_n^2
&= \Big( \frac{b-a}{2} \Big)^2 + \frac{1}{n} \cdot \Big( \frac{b-a}{2} \Big)^2 \quad \quad \quad \quad \quad \ \\[6pt]
&= \Big( 1 + \frac{1}{n} \Big) \cdot \Big( \frac{b-a}{2} \Big)^2 \\[6pt]
&= \frac{n+1}{n} \cdot \Big( \frac{b-a}{2} \Big)^2 \\[6pt]
&= \frac{(n-1)(n+1)}{n^2} \cdot \frac{n}{n-1} \cdot \Big( \frac{b-a}{2} \Big)^2 \\[6pt]
&= \frac{n^2-1}{n^2} \cdot \frac{n}{n-1} \cdot \Big( \frac{b-a}{2} \Big)^2. \\[6pt]
\end{align}$$
This is a slightly lower sample variance value than when we have an even number of data points, but the two cases converge when $n \rightarrow \infty$.