Denote the support by $[a,b]$ and the variance by $\sigma^2$. Repeating the variational calculus argument that shows the normal distribution maximizes the differential entropy (see here for instance), probability distributions of variance $\sigma^2$ and supported in $[a,b]$ that attain the largest possible entropy should have a density function of the form

$$
p(x)=\frac{e^{-\lambda(x-\mu)^2}}{\int_a^be^{-\lambda(t-\mu)^2}{\rm{d}}t}\quad
(a\leq x\leq b)\tag{$\star$}
$$
where $\lambda$ and $\mu$ satisfy

$$
\begin{cases}
\int_a^bte^{-\lambda(t-\mu)^2}{\rm{d}}t=\mu\int_a^be^{-\lambda(t-\mu)^2}{\rm{d}}t,\\
\int_a^b(t-\mu)^2e^{-\lambda(t-\mu)^2}{\rm{d}}t=\sigma^2\int_a^be^{-\lambda(t-\mu)^2}{\rm{d}}t.
\end{cases}
$$
The first equation holds automatically if $\mu$ is at the center of the interval $(a,b)$. Assuming that $a,b$ are finite and the variational problem has a unique solution, then $\mu$ must be $\frac{a+b}{2}$: If $p:[a,b]\rightarrow(0,\infty)$ is the density function for the maximizer, then $x\in [a,b]\mapsto p(a+b-x)$ defines another probability distribution with the same support, variance and differential entropy. Then, the first equation above becomes irrelevant and one only needs to solve
$$
\int_a^b\left(t-\frac{a+b}{2}\right)^2e^{-\lambda\left(t-\frac{a+b}{2}\right)^2}{\rm{d}}t=\sigma^2\int_a^be^{-\lambda\left(t-\frac{a+b}{2}\right)^2}{\rm{d}}t \tag{$\star\star$}
$$
for $\lambda$.

**Update)** Following the comment by KConrad, a more general system of equations should be considered. Let's write the density function as

$$
p(x)=\frac{e^{-(\lambda x^2+\beta x)}}{\int_a^be^{-(\lambda x^2+\beta x)}{\rm{d}}t}\quad
(a\leq x\leq b), \tag{$\star\star\star$}
$$
In $(\star)$, it was implicitly assumed that $-\frac{\beta}{2\lambda}$ coincides with the mean $\mu$. Given the form $(\star\star\star)$ for the density function, the system of equations becomes

$$
\begin{cases}
\int_a^bte^{-(\lambda t^2+\beta t)}{\rm{d}}t=
\mu\int_a^be^{-(\lambda t^2+\beta t)}{\rm{d}}t,\\
\int_a^b(t-\mu)^2e^{-(\lambda t^2+\beta t)}{\rm{d}}t=
\sigma^2\int_a^be^{-(\lambda t^2+\beta t)}{\rm{d}}t.
\end{cases}
$$

As argued before, assuming that the solution is unique, $p(x)$ must coincide with $p(a+b-x)$. When $p(x)$ is of the form $(\star\star\star)$, by comparing the coefficients of $x$ in the exponent from the numerator, $p(a+b-x)\equiv p(x)$ implies $-\frac{\beta}{2\lambda}=\frac{a+b}{2}$. Moreover,
$\int_a^btp(t){\rm{d}}t=\int_a^btp(a+b-t){\rm{d}}t$ implies that the mean $\mu=\int_a^btp(t){\rm{d}}t$ is equal to $\frac{a+b}{2}$. Just like before, when $\mu=-\frac{\beta}{2\lambda}=\frac{a+b}{2}$, the first equation holds automatically because
$\int_a^b\left(t-\frac{a+b}{2}\right)e^{-\lambda\left(t-\frac{a+b}{2}\right)^2}{\rm{d}}t=0$. So it only remains the second equation which would be same as $(\star\star)$, and should be solved for $\lambda$ to derive a density function of the form
$$
p(x)=\frac{e^{-\lambda\left(x-\frac{a+b}{2}\right)^2}}{\int_a^be^{-\lambda\left(t-\frac{a+b}{2}\right)^2}{\rm{d}}t}.
$$