The purpose of the present paper is to incorporate stochastic interest rates into a matrix-approach to multi-state life insurance, where formulas for reserves, moments of future payments and equivalence premiums can be obtained as explicit formulas in terms of product integrals or matrix exponentials. To this end we consider the Markovian interest model, where the rates are piecewise deterministic (or even constant) in the different states of a Markov jump process, and which is shown to integrate naturally into the matrix framework. The discounting factor then becomes the price of a zero-coupon bond which may or may not be correlated with the biometric insurance process. Another nice feature about the Markovian interest model is that the price of the bond coincides with the survival function of a phase-type distributed random variable. This, in particular, allows for calibrating the Markovian interest rate models using a maximum likelihood approach to observed data (prices) or to theoretical models like e.g. a Vasiček model. Due to the denseness of phase-type distributions, we can approximate the price behaviour of any zero-coupon bond with interest rates bounded from below by choosing the number of possible interest rate values sufficiently large. For observed data models with few data points, lower dimensions will usually suffice, while for theoretical models the dimensionality is only a computational issue.