Experience working with distributions suggests analyzing the logarithm of the density function, rather than the density itself. Because the log is a monotonic increasing transformation, the mode of the log density occurs at the same value as the mode of the density. (This approach has general application, not just for beta distributions.)
Let's develop this solution in steps, asking Mathematica to simplify as much as possible at each stage. First, the distribution itself:
f[x_, {a_, b_}] := PDF[BetaDistribution[a, b]][x]
(Bear in mind that this is defined only for positive values of $a$ and $b$.)
Now, its logarithm. Take care to avoid places where the distribution will be zero or undefined:
g[x_, {a_, b_}] := Evaluate[Simplify[Log[f[x, {a, b}]], Assumptions -> 0 < x < 1]]
We plan to inspect the zeros of its derivative for the mode:
h[x_, {a_, b_}] := Evaluate[Simplify[D[g[x, {a, b}], {x}]]]
The critical points will include all zeros of the derivative:
criticalPoints = Solve[h[x, {a, b}] == 0, {x}]
$\left\{\left\{x\to \frac{-1+a}{-2+a+b}\right\}\right\}$
If we're working blindly with an unfamiliar distribution, we had better check the endpoints to make sure that one of them isn't actually a mode, too. This can be problematic, because Mathematica often will have trouble finding the limiting values for symbolic parameters $a$ and $b$. One alternative is to explore the situation graphically and manipulate $a$ and $b$ over reasonable ranges:
Manipulate[
TableForm[{"Value at 0: " <> ToString[Limit[f[x, {a, b}], x -> 0, Direction -> -1]],
"Value at 1: " <> ToString[Limit[f[x, {a, b}], x -> 1, Direction -> 1]],
Plot[f[x, {a, b}], {x, 0, 1}, ImageSize -> 400]}],
{{a, 1/2}, 0, 2}, {{b, 3/2}, 0, 2}]

This will quickly show that
When $a \lt 1$, there is a mode at $0$.
When $b \lt 1$, there is a mode at $1$.
Otherwise, when both $a \ge 1$ and $b \ge 1$, the mode is at $\frac{a-1}{[a-1] + [b-1]}$, as given by the value of criticalPoints
.
When both $a \lt 1$ and $b \lt 1$, the distribution is bimodal ("U shaped").
Commonest
gives the mode of a list, so you could do this on a sample, but I assume what you really want is to get the maximum of thePDF
by taking the derivative and solving for zero? $\endgroup$mode=ArgMax[PDF[your distribution, x],x]
: for example,ArgMax[PDF[BetaDistribution[2, 3], x], x]
gives1/3
. $\endgroup$