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 → ∞, 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 = 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 × 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.