There are indeed $\binom{52}{13}$ different 13-card hands and this will indeed be the size of our sample space and thus our denominator when we finish our calculations.
For the numerator, we need to pause for a moment and understand what the problem is actually asking, since this appears to be where you got stuck.
We are asked to find the probability that in our hand of thirteen cards, there is at least one suit for which we have all three face cards. For example $(A\spadesuit,2\spadesuit,3\spadesuit,\dots,10\spadesuit,J\spadesuit,Q\spadesuit,K\spadesuit)$ has all three of the face cards for spades. Similarly if all those cards happened to be hearts instead it would also count since we would have all of the face cards for hearts. Similarly still, a hand like $(J\spadesuit,Q\spadesuit,K\spadesuit,J\heartsuit,Q\heartsuit,K\heartsuit,J\diamondsuit,Q\diamondsuit,K\diamondsuit,\dots)$ would count since we have all of the face cards from spades (we also happen to have all of the face cards from hearts and diamonds too).
Let $\spadesuit$ represent the event that we have have all of the face cards from spades. Similarly, let $\diamondsuit, \heartsuit, \clubsuit$ represent the event that we have all of the face cards from diamonds, hearts, and clubs respectively.
You are asked to find $Pr(\spadesuit\cup\diamondsuit\cup \clubsuit\cup \heartsuit)$
To do this, let us apply inclusion exclusion. We expand the above as:
$Pr(\spadesuit\cup \diamondsuit\cup\clubsuit\cup\heartsuit) = Pr(\spadesuit)+Pr(\diamondsuit)+\dots-Pr(\spadesuit\cap \diamondsuit)-Pr(\spadesuit\cap \clubsuit)-\dots+Pr(\spadesuit\cap \diamondsuit\cap \clubsuit)+\dots-Pr(\spadesuit\cap\diamondsuit\cap \clubsuit\cap \heartsuit)$
Now, let us calculate each individual term in the expansion.
The calculation you did before is relevant. Indeed, we calculate $Pr(\spadesuit)=\dfrac{\binom{3}{3}\binom{49}{10}}{\binom{52}{13}}$. This is again merely the probability that we have all of the face cards from the spades and is not the final probability that we were tasked with calculating.
We continue and calculate more terms:
For example $Pr(\spadesuit\cap \diamondsuit)=\dfrac{\binom{6}{6}\binom{46}{7}}{\binom{52}{13}}$
We then notice what symmetry there is in the terms and can simplify some. Finally, we write the final expression for our final answer (and get an exact number only if actually requested or required, usually opting to leave the answer in terms of binomial coefficients without additional simplification).
A simple estimate: given $n$ random variables, independently and identically distributed evenly over $M$ states, the probability that some two of them are equal is at most $\frac{n(n-1)}{2M}$. Why? Because that's the expected value for the number of pairs that are equal; there are $\binom{n}{2}$ pairs, each with a probability of $\frac1M$ of being equal.
So then, if $n$ is significantly less than $\sqrt{M}$, the probability of some two being equal is small.
Now, we focus on the shuffling problem. How big is $52!$, really? For that, we have Stirling's formula:
$$n!\approx \frac{n^n}{e^n}\cdot\sqrt{2\pi n}$$
Take logarithms, and we get $n(\log n-\log e)+\frac12\log(2\pi n)$. Calculating that in the spreadsheet I've got open already, I get a base 10 log of about $67.9$, for about $8\cdot 10^{67}$ possible shuffled decks.
So, then, if we take $10^{30}$ shuffles, that's a probability of less than $\frac{5\cdot 10^{59}}{8\cdot 10^{67}}< 10^{-8}$ (one in a hundred million) that any two of them match exactly. That $10^{30}$ - if you had a computer for every person on earth, each generating shuffled decks a billion times per second, the planet would be swallowed up by the sun before you got there.
Best Answer
Your original answer of $\dfrac{3 \times 10^{14}}{52!}$ is not far from being right. That is in fact the expected number of times any ordering of the cards has occurred.
The probability that any particular ordering of the cards has not occurred, given your initial assumptions, is $\left(1-\frac1{52!}\right)^{(3\times10^{14})}$, and the probability that it has occurred is 1 minus this value. But for small values of $n\epsilon$, $(1+\epsilon)^n$ is nearly $1+n\epsilon$. In particular, since $52!\approx 8\times 10^{67}$ and so $\dfrac{3\times10^{14}}{52!}\approx 3.75\times 10^{-54}$ is microscopically small, $1-\left(1-\frac1{52!}\right)^{(3\times10^{14})}$ is very nearly $\frac1{52!}\times (3\times10^{14})$.