As already noted in the comments, a possible initial approach to this problem is the following. Let us suppose that, after selecting the maximal number in each row, the minimum number $A$ is that in the $j^{th} $ row. Also, let us suppose that, after selecting the minimum number in each column, the maximal number $B$ is that in the $i^{th} $ column. Now consider the number $x_{i,j} $ corresponding to the crossing point of the $i^{th} $ column and $j^{th} $ row. We directly get that $A \geq x_{i,j} \geq B $. Thus, the searched probability that $A >B $ is equal to $1-Pr[ A = x_{i,j} = B] $, where this second term expresses the probability that both the two procedures finally identify the same number in the matrix.
We can now continue as follows. First, note that the condition that both two procedures finally identify the same number in the matrix implies that there exists a number $x_{i,j} $ in the matrix which is the highest in its row, and the lowest in its column. Also note that, if such a number exists, then it must be unique. To show this, let us assume that there exists another number $x_{r,s} $ (corresponding to the crossing point of the $r^{th} $ column and $s^{th} $ row) that - just like $x_{i,j} $ - is the highest in its row, and the lowest in its column. This would imply $$x_{i,j} > x_{r,j} > x_{r,s} > x_{r,j} > x_{i,j} $$ which is clearly impossible. Therefore, we have that $Pr[ A = x_{i,j} = B]$ is equivalent to the probability that there exists a number $x_{i,j} $ in the matrix which is the highest in its row, and the lowest in its column.
We can now try to calculate this probability. The probability $P_C(k)$ that, after dividing the first $n^2$ integers in $n $ random groups of $n $ elements (i.e., the columns), a given number $k$ of this set is the lowest in its group is given by
$$P_C(k)= \frac {n^2-k}{n^2-1} \cdot \frac{n^2-k-1}{n^2-2}.... \cdot \frac {n^2-k-n+2} {n^2-n+1} $$
where the sequence of fractions expresses the probability that, given $k $, then the first, the second... and the $(n-1)^{th}$ among the other numbers in its group/column are all $>k $. This probability formula is valid only for $k \leq n^2-n+1$ (for higher values of $k $ the probability is zero). The expression above can also be written as
$$ P_C(k)= \frac {(n^2-k)!}{(n^2-k-n+1)!} \cdot \frac {(n^2-1)!}{ (n^2-n)!} =\frac {\binom{n^2-k}{n-1}} {\binom{n^2-1}{ n-1}} $$
By similar considerations, we can calculate the probability $ P_R(k) $ that, after dividing the first $n^2$ integers in other $n $ random groups of $n $ elements (i.e., the rows), our given number $k$ is the highest in its group. Because this time we have to exclude, from the possible other terms that can compare in this group, the $n-1$ terms already considered in the same column of $k $ (it is clear that these terms cannot compare also in the same row of $k $), and taking into account that all these $n-1$ terms are $>k $, we obtain
$$P_R(k) = \frac {k-1}{n^2-n} \cdot \frac{k-2}{n^2-n-1}.... \cdot \frac {k-n+1} {n^2-2n+2} $$
where again the sequence of fractions expresses the probability that, given $k $, then the first, the second... and the $(n-1)^{th}$ among the other numbers in its group/column are all $<k $. This probability formula is valid only for $k \geq n$ (for lower values of $k $ the probability is zero). Note that, for the reason explained above, the denominators have been decreased by $n-1$ as compared to those used before. The expression above can be written as
$$ P_R(k) = \frac {(k-1)!}{(k-n)!} \cdot \frac {(n^2-n)!}{ (n^2-2n+1)!} =\frac {\binom{k-1}{ n-1}} {\binom{n^2-n}{ n-1}} $$
We can now calculate the probability that in the matrix there exists a number $x_{i,j} $ that is the highest in its row and the lowest in its column. This can be obtained by summing, for all possible values of $k $, the product of $P_C(k) $ and $P_R(k)$:
$$\sum _{k=n}^{n^2-n+1} P_C(k) \cdot P_R(k) = \sum _{k=n}^{n^2-n+1} \frac {\binom{n^2-k}{ n-1}} {\binom{n^2-1}{ n-1} } \frac {\binom {k-1}{ n-1}} {\binom{n^2-n}{ n-1} } $$
Reminding that this probability is equal to $Pr[ A = x_{i,j} =B]$, we finally obtain
$$Pr[ A>B] = 1- Pr[A= x_{i,j}=B] $$
$$=1- \sum _{k=n}^{n^2-n+1} \frac {\binom{n^2-k}{ n-1}} {\binom{n^2-1}{ n-1} } \frac {\binom{k-1}{ n-1}} {\binom{n^2-n}{ n-1}} $$
Note that, for the case $n=2$, this expression reduces to the cases $k=2$ and $k=3$. Since the summation gives $1/3$ in both cases, the final result for the case $n=2$ is $$Pr [A>B]=1-1/3-1/3=1/3$$ as anticipated in the comments. The cases $n=4$ and $n=5$ give $$1-3/10=7/10$$ and $$1-4/35=31/35$$
respectively, which are very near to the experimental values reported in the OP. The probability rapidly grows up and tends to $1$. For example, for $n=10$ its value is $$1-5/46189 \approx 0.99989... $$
again confirming the experimental value.
To get a closed value, we can transform the binomial coefficients using factorials and simplify, so that the final result becomes
$$Pr [A >B] =1- \frac {2n \, (n!)^2}{(2n!)}$$
The rapid increase of the function and its first values for small $n $ are shown by WA here.
