The random variable 1 + z1 + z1z2 +. .. appears in many contexts and was shown by Kesten to exhibit a heavy tail distribution. We consider natural extensions of this variable and its associated recursion to N × N matrices either real symmetric β = 1 or complex Hermitian β = 2. In the continuum limit of this recursion, we show that the matrix distribution converges to the inverse-Wishart ensemble of random matrices. The full dynamics is solved using a mapping to N fermions in a Morse potential, which are non-interacting for β = 2. At finite N the distribution of eigenvalues exhibits heavy tails, generalizing Kesten's results in the scalar case. The density of fermions in this potential is studied for large N , and the power-law tail of the eigenvalue distribution is related to the properties of the so-called determinantal Bessel process which describes the hard edge universality of random matrices. For the discrete matrix recursion, using free probability in the large N limit, we obtain a self-consistent equation for the stationary distribution. The relation of our results to recent works of Rider and Valkó, Grabsch and Texier, as well as Ossipov, is discussed.