What are the properties of the distribution of medians from samples of a set of integers? Given a set of integers $S = \{1, 2, 3, ..., n\}$ I draw uniformly at random without repetition $m$ integers a large number of times ($1 \leq m \leq n$). For each of these samples (each of size $m$) I compute the median. What is the expected distribution of these medians?
By intuition, we can say the sample space from which we can take the samples is the set of all permutations without repetition of the elements of $S$ of size $m$. Different permutations of this set will have the same median making certain median values more likely than others so the resulting distribution is not uniform for $1 \lt m \lt n$. For the extreme cases $m = 1$ we should observe a uniform distribution in the range $[1, n]$ and for $m = n$ all medians will be the median of $S$. I am unsure on what to state about the distribution of medians for all other values though.

These order statistics follow a shifted beta-binomial distribution
Let $X_1,...,X_m \sim \text{SRSWOR}(n)$ denote the sampled values.  If you have an odd number of values in the sample then the median is just the middle order statistic.  I am going to generalise from this and show how you can obtain the distribution of any order statistic $X_{(k)}$ over the order values $k =1,...,m$.  (Getting the median in the case where you have an even number of values is a bit trickier, since it involves the joint distribution of two adjacent order statistics; I will leave this as an exercise.)
There are a couple of different ways you can derive the distribution for this order statistic.  The simplest method is through a direct combinatorial argument.  There are ${n \choose m}$ possible samples that can be drawn from the population (with equal probability).  To get the event $X_{(k)} = x$ we have ${x-1 \choose k-1}$ ways of choosing the previous order statistics and ${n-x \choose m-k}$ ways of choosing the subsequent order statistics.  Consequently, we get a shifted beta-binomial distribution:
$$\mathbb{P}(X_{(k)} = x)
= \frac{{x-1 \choose k-1} {n-x \choose m-k}}{{n \choose m}}
= \text{BetaBin}(x-k|n-m,k,m-k+1).$$
(To confirm that this is a valid probability mass function ---i.e., that it sums to one over the support--- note this related summation result for products of binomial coefficients.)  It is also possible to obtain this result by using the relationship of the order statistics to the hypergeometric distribution, but this latter method is much more complicated (omitted here).
With a bit of algebra, it is possible to derive the general form for the factorial moments of this distribution, and through this obtain the mean and variance.  These are given by:
$$\mathbb{E}(X_{(k)}) = \frac{n+1}{m+1} \cdot k
\quad \quad \quad \quad \quad 
\mathbb{V}(X_{(k)}) = \frac{(n+1)(n-m)}{(m+1)^2 (m+2)} \cdot k(m-k+1).$$
Further examination of this distribution can be found in O'Neill (2022).  The distribution of the scaled order statistics can be shown to have a mixture distribution that is closely related to the order statistics from IID observations from a continuous uniform distribuion on the unit interval.

Implementation in R: The beta-binomial distribution is available in various packages in R (e.g., the extraDistr and TailRank packages).  It is also quite simple to program the specific variant shown above (using the required parameterisation and shift) to give the distribution of the order statistic, which I will give below as the FPOS distribution (standing for finite-population-order-statistic).  Appealing to standard notation in these problems, I will use N for the population size and n for the sample size.  The density function is programmed using computation in log-space, in order to ensure computational stability.
#Probability mass function for FPOS distribution
dFPOS <- function(x, k, n, N, log = FALSE) {

  #Check input types
  if (!is.numeric(x))    stop('Input x must be numeric')
  if (!is.numeric(k))    stop('Order k must be a positive integer')
  if (!is.numeric(n))    stop('Sample size n must be a positive integer')
  if (!is.numeric(N))    stop('Population size N must be a positive integer')
  if (!is.logical(log))  stop('Input log must be a logical value')
  ki <- as.integer(k)
  ni <- as.integer(n)
  Ni <- as.integer(N)
  if (k != ki)           stop('Order k must be a positive integer')
  if (n != ni)           stop('Sample size n must be a positive integer')
  if (N != Ni)           stop('Population size N must be a positive integer')

  #Check input lengths and values
  if (length(k) != 1)    stop('Order k should be a single positive integer')
  if (length(n) != 1)    stop('Sample size n should be a single positive integer')
  if (length(N) != 1)    stop('Population size N should be a single positive integer')
  if (length(log) != 1)  stop('Input log should be a single logical value')
  if (k < 1)             stop('Order k must be positive')
  if (n < 1)             stop('Sample size n must be positive')
  if (N < 1)             stop('Population size N must be positive')
  if (k > n)             stop('Order k cannot be greater than sample size n')
  if (n > N)             stop('Sample size n cannot be greater than population size N')
  
  #Compute log-probabilities
  LOGPROBS <- rep(0, N)
  for (t in 1:N) {
    LOGPROBS[t] <- lchoose(t-1, k-1) + lchoose(N-t, n-k) - lchoose(N, n) }
  
  #Generate output
  L <- length(x)
  OUT <- rep(-Inf, L)
  for (i in 1:L) {
    xx <- x[i]
    if (xx %in% (1:N)) {
      OUT[i] <- LOGPROBS[xx] } }
  
  #Return output
  if (log) { OUT } else { exp(OUT) } }

Using this function we can obtain the distribution of the median (for an odd sample size) by taking $m=2k-1$ and finding the distribution of the order statistic $X_{(k)}$.  In the code below we compute and plot the mass function for the median for the case where $n=16$ and $m=9$.
#Set parameters
n <- 16
m <- 9

#Compute and plot the mass function
PROBS <- dFPOS(1:n, k=(m+1)/2, N = n, n = m)
names(PROBS) <- 1:n
barplot(PROBS, xlab = 'Value', ylab = 'Probability', col = 'blue')


