A multi-state life insurance model is described naturally in terms of the intensity matrix of an underlying (time-inhomogeneous) Markov process which specifies the dynamics for the states of an insured person. Between and at transitions, benefits and premiums are paid, defining a payment process, and the technical reserve is defined as the present value of all future payments of the contract. Classical methods for finding the reserve and higher order moments involve the solution of certain differential equations (Thiele and Hattendorff, respectively). In this paper we present an alternative matrix-oriented approach based on general reward considerations for Markov jump processes. The matrix approach provides a general framework for effortlessly setting up general and even complex multi-state models, where moments of all orders are then expressed explicitly in terms of so-called product integrals of certain matrices. Thiele and Hattendorff type of theorems may be retrieved immediately from the matrix formulae. As a main application, methods for obtaining distributions and related properties of interest (e.g. quantiles or survival functions) of the future payments are presented from both a theoretical and practical point of view, employing Laplace transforms and methods involving orthogonal polynomials.