  We study a model of $n$ non-intersecting squared Bessel processes in the
confluent case: all paths start at time $t = 0$ at the same positive value $x =
a$, remain positive, and are conditioned to end at time $t = T$ at $x = 0$. In
the limit $n \to \infty$, after appropriate rescaling, the paths fill out a
region in the $tx$-plane that we describe explicitly. In particular, the paths
initially stay away from the hard edge at $x = 0$, but at a certain critical
time $t^*$ the smallest paths hit the hard edge and from then on are stuck to
it. For $t \neq t^*$ we obtain the usual scaling limits from random matrix
theory, namely the sine, Airy, and Bessel kernels. A key fact is that the
positions of the paths at any time $t$ constitute a multiple orthogonal
polynomial ensemble, corresponding to a system of two modified Bessel-type
weights. As a consequence, there is a $3 \times 3$ matrix valued
Riemann-Hilbert problem characterizing this model, that we analyze in the large
$n$ limit using the Deift-Zhou steepest descent method. There are some novel
ingredients in the Riemann-Hilbert analysis that are of independent interest.
