I don't know if this counts as an elegant solution in your book, but I think it's cute.
Let's say the "frequency state" of a deck is the number of cards of each face value remaining. A full deck, for example, has the frequency state "4 aces, 4 twos, 4 threes...," while an empty deck has the frequency state "0 aces, 0 twos, 0 threes...." There are $5^{13}$ possible frequency states. When you draw a card from a deck, the frequency state changes in a way that depends only on the face value of the card you drew.
You can turn the set of possible frequency states into a directed graph by drawing an arrow for each way the frequency state can change when you draw a card. I'll call this graph the "draw graph." Each vertex in the draw graph has at most 13 edges leaving it, one for each type of card you could draw.
You can turn the draw graph into a weighted directed graph by assuming that cards are drawn uniformly at random, and weighting each arrow by the probability of that kind of draw. The full-deck state, for example, has 13 arrows leaving it, each with weight $\tfrac{1}{13}$. If you draw a queen, you end up in a state that still has 13 arrows leaving it, but 12 of them have weight $\tfrac{4}{51}$, and one—the arrow for drawing another queen—has weight $\tfrac{3}{51}$. The weights of the arrows leaving each state add up to one, so the draw graph is a Markov chain.
Let's say the draw has been passed to the dealer. We know the dealer's hand and the frequency state of the deck. Here's a cool fact: from now on, we can figure out the dealer's hand just by looking at the frequency state of the deck. That's because, when the dealer starts hitting, all the cards she draws from the deck end up in her hand. Using this fact, we can translate properties of the dealer's hand, like whether it's bust, into properties of the frequency state.
Let's record this information on the draw graph by labeling its vetrtices. We'll label the states where the dealer stays as "stay states," the states where the dealer is bust as "bust states," and the states where the dealer keeps hitting as "hit states." When the dealer is hitting, she's walking randomly along the draw graph, with her direction from each state chosen using the arrow weights as probabilities. She keeps walking until she reaches a stay state or a bust state.
The dealer has to eventually stay or go bust, so the process we just described is an absorbing Markov chain. Like most things in graph theory, absorbing Markov chains have a very pretty linear algebraic description in terms of the adjacency map of the transition graph. If you know how this works, I can describe very quickly how to calculate the bust probability.
Cribbing from Wikipedia, let $Q$ be the map describing transitions from hit states to hit states, and let $R$ be the map describing transitions from hit states to stay and bust states. Let $\pi$ be the vector corresponding to the initial state, and let $\Pi$ be the projection map onto the subspace spanned by the bust states. The bust probability is the sum of the entries of the vector
$$\Pi R (1 - Q)^{-1}\pi.$$
The $(1 - Q)^{-1}$ factor describes the part of the process where the dealer hits over and over, so I think it encapsulates the tricky "recursive bit" of the computation.
(Caution: my operators are the transposes of Wikipedia's, which is why the formula looks backwards.)
I think this question is related to a broader question that I ask myself all the time in research: what does it mean to have a "nice solution" to a problem? When I was young, I was taught that a "nice solution" is a formula for the thing you want to calculate, but that's not always true! Having a forumla for something often tells you very little about it, and other descriptions are often much more useful from a practical point of view.
I'm not sure whether the description of the bust probability given above is much use, but for this problem, I suspect that a linear algebraic description of this kind will be more useful than a formula.
Best Answer
This is my second try, but I think is a MUCH better argument than what I had before. So I deleted the more complicated answer I tried to give earlier.
The total # of possible deals to the two players is: $$\binom{52}{2}\binom{50}{2}=1,624,350$$ (i.e., choose two cards for the first player and then two for the second.)
Then, the # of ways to deal a blackjack to both players would be: $$\binom{4}{1}\binom{16}{1}\binom{3}{1}\binom{15}{1}=2880$$ (i.e, choose which ace to give to the first player, then which 10,J,Q,K for the first player, then which ace for the second player, and then which 10,J,Q,K for the second player.)
Also, the # of ways that exactly one player gets a blackjack would be: $$\binom{2}{1}\binom{4}{1}\binom{16}{1}\times(\binom{3}{2}+\binom{3}{1}\binom{32}{1}+\binom{47}{2})=151,040$$ (i.e., choose which player to give the blackjack, then choose which ace to give them, and then which 10,J,Q,K to give them. Then for the other player we either give them two aces from the 3 that are left, give them one ace from the 3 that are left and 1 of the 32 2-9 cards, or give them two non-aces from the 47 non-aces that are left.)
So, the probability of at least one player getting a blackjack is: $$\frac{2880+151,040}{1,624,350}$$
Hence the probability of neither player getting a blackjack is: $$1- \frac{2880+151,040}{1,624,350} \approx 90.5\%$$