Your calculation is correct, and you were actually quite close to the more elegant answer yourself.
You wrote:
"It seems that simply $n\cdot0.00003232$ will give a good
approximation to the probability I'm after but something seems off
about this in that $1$ player receiving a royal flush would negate
anyone else from receiving a royal flush (except in the case of a
royal being on the board)."
The opposite is actually true: If one player receiving a royal flush would prevent anyone else from receiving one, the events would be disjoint and the probability for any player to get a royal flush would indeed be simply $n$ times the probability for a particular player to get one. As you noticed, there is only a single case in which more than one player has a royal flush, namely if it's on the board, and we have to compensate for overcounting this case. In adding the $n$ probabilities as if they were disjoint, we're counting this case $n$ times instead of once, so we need to subtract it $n-1$ times. A royal flush being on the board is one of $\binom75=21$ ways of having a royal flush in $7$ cards. Thus the probability including the correction for overcounting is
$$
np-\frac1{21}(n-1)p=\frac{20n+1}{21}p\;,
$$
where
$$
p=\frac{4\binom{47}2}{\binom{52}7}=\frac1{30940}
$$
is the probability that you calculated for a $7$-card hand to have a royal flush. Thus the probability you want is
$$
\frac{20n+1}{649740}\;,
$$
in agreement with your calculation.
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).
Best Answer
A start: Line up the players in order of student number, and call them Player 1, Player 2, and so on up to Player 13.
There are $\binom{52}{4}$ ways to choose the cards that Player 1 gets, and for each of these ways there are $\binom{48}{4}$ ways to choose the cards that Player 2 gets, and so on, for a total of $\binom{52}{4}\binom{48}{4}\binom{44}{4}\cdots \binom{4}{4}$ ways.
This number simplifies to $\dfrac{52!}{(4!)^{13}}$. All these ways are equally likely.
For event $A$, we count the "favourables." The spades that the various players get can be chosen in $13!$ ways. For each of these ways, the hearts they get can be chosen in $13!$ ways, and so on, for a total of $(13!)^4$.
Next we count the "favourables" for event $B$. The kinds the various players get can be chosen in $13!$ ways.
Remark: The exactly $3$ will be tricky. It is not difficult to give an expression for the number of ways to distribute the "threes." However, for determining the distributions of the "odd" cards, we will need to ensure that everybody's odd card is of a different kind than the main $3$. To count, we need to know something about derangements.