In view of the [algorithms] tag (and since you are a software engineer), perhaps you'll be satisfied with the following answer. Assume that $X_i$ are independent ${\rm Beta}(\alpha_i,\beta_i)$ variables. Then, you can evaluate the probability ${\rm P}(X_i = \max \lbrace X_1 , \ldots ,X_n \rbrace )$ using Monte Carlo simulations, as follows. Obviously, the problem amounts to simulating a ${\rm Beta}(\alpha,\beta)$ variable. This can be done simply as follows, according to Example 2.11 in the book Monte Carlo statistical methods (see references therein). If $U$ and $V$ are independent ${\rm uniform}[0,1]$ variables, then the distribution of
$\frac{{U^{1/\alpha } }}{{U^{1/\alpha } + V^{1/\beta } }}$ conditional on $U^{1/\alpha } + V^{1/\beta } \le 1$ is the ${\rm Beta}(\alpha,\beta)$ distribution. As noted in that example, this result does not provide a good algorithm to generate ${\rm Beta}(\alpha,\beta)$ variables for large values of $\alpha$ and $\beta$ (because of the constraint on $U^{1/\alpha } + V^{1/\beta }$). But if your $\alpha_i$ and $\beta_i$ are not large, you might find this simple algorithm useful enough (depending on the accuracy you wish to achieve).

**EDIT**: This approach may be particularly useful for values $\alpha_i,\beta_i \in (0,1)$, for two reasons. First, this increases the probability that $U^{1/\alpha } + V^{1/\beta } \le 1$ (that is, the pair $(U,V)$ is not rejected). Second, the ${\rm Beta}(\alpha,\beta)$ density function is not bounded if $\alpha \in (0,1)$ or $\beta \in (0,1)$, and so a tractable analytical expression for ${\rm P}(X_i = \max \lbrace X_1 , \ldots ,X_n \rbrace )$ is not likely to be found in this case. Of course, everything changes if the parameters are integers...