Can Mathematica find an expression for the distribution of the median of N i.i.d. random variables?

Mathematica does make this pretty easy. The statistic of interest is the typical estimator of the median when the sample size is even. When the sample size is odd the sample median has a beta distribution:

OrderDistribution[{UniformDistribution[{0, 1}], n}, (n + 1)/2]
(* BetaDistribution[(1 + n)/2, 1 + 1/2 (-1 - n) + n] *)

Now for the case when $n$ is even. First find the joint distribution of the middle two order statistics. Then find the distribution of the mean of those two statistics.

n = 6;
od = OrderDistribution[{UniformDistribution[{0, 1}], n}, {n/2, n/2 + 1}];

md = TransformedDistribution[(x1 + x2)/2, {x1, x2} \[Distributed] od];

PDF[md, x]

Plot[Evaluate[PDF[md, x]], {x, 0, 1}]

To obtain the distribution for general $n$ when $n$ is even we have to use some other than TransformedDistribution. We need to integrate the joint density function and treat $0<x<1/2$, $x=1/2$, and $1/2<x<1$ separately.

fltOneHalf = 2 Integrate[(x1^(-1 + n/2) (1 - x2)^(-1 + n/2) n!)/((-1 + n/2)!)^2 /. 
  x2 -> 2 x - x1, {x1, 0, x}, Assumptions -> n > 1 && 0 < x < 1/2]
(* -((4 ((1 - 2 x) x)^(n/2) Gamma[n]*
  Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, x/(-1 + 2 x)])/((-1 + 2 x)*
  Gamma[n/2]^2)) *) 

fOneHalf = 2 Integrate[(x1^(-1 + n/2) (1 - x2)^(-1 + n/2) n!)/((-1 + n/2)!)^2 /.
  x2 -> 1 - x1, {x1, 0, 1/2}, Assumptions -> n > 1]
(* (2^(2 - n) n!)/((-1 + n) ((-1 + n/2)!)^2) *)

(* Because the density is symmetric, we'll take advantage of that *)
fgtOneHalf = FullSimplify[fltOneHalf /. x -> y /. y -> 1 - x]
(* (4 (-1 + (3 - 2 x) x)^(n/2) Gamma[n]*
  Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, (-1 + x)/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2) *)

Putting this together in a single function:

pdf[n_, x_] := 
 Piecewise[{{-((4 ((1 - 2 x) x)^(n/2)*Gamma[n] Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, 
       x/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2)), 0 < x < 1/2},
   {(2^(2 - n) n!)/((-1 + n) ((-1 + n/2)!)^2), x == 1/2},
   {(4 (-1 + (3 - 2 x) x)^(n/2) * Gamma[n]*
 Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, (-1 + x)/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2), 
    1/2 < x < 1}}, 0]

Here's an interesting counter-example for a discrete uniform distribution does not tend to your shape as $N$ grows.

Let your r.v. $x$ be distributed as per a coin toss, taking value ${0,1}$ with equal probability. Then you have three possible outcomes after $N$ coin tosses; $0, 1/2$ , or $1$

The middle outcome's probability is the probability that with $N$ coin tosses you get exactly $N/2$ zeros or ones. This probability is

$$ 2^{-n} \binom{n}{\frac{n}{2}} $$

Which decreases in $N$, resulting in the plot below, with the distribution of probabilities of values with $N$

EDIT

You can see the same effect with the discrete distribution for $x$ expanded to take values $x=\{1, 2, ... , 50\}$ with equal probability. The non-integer values are much less likely, since the odds of the two middle points hitting the boundary is low. The integers values do tend to a Guassian.

middleMean[n_, range_] := Module[{res}, 
     res = Sort@Table[RandomChoice[Range[range]], {n}];
     Mean[{res[[n/2]], res[[n/2 + 1]]}]
  ]

Histogram[Table[middleMean[500, 50], {1000}], 50]