A common concern with Markov reliability and reward models is that model parameters, i.e., component failure and repair rates, are seldom perfectly known. This paper proposes a numerical method based on the Taylor series expansion of the underlying Markov chain stationary distribution (associated to the reliability and reward models) to propagate parametric uncertainty to reliability and performability indices of interest. The Taylor series coefficients are expressed in closed form as functions of the Markov chain generator-matrix group inverse. Then, to compute the probability density functions of the reliability and performability indices, random variable transformations are applied to the polynomial approximations that result from the Taylor series expansion. Additionally, closed-form expressions that approximate the expectation and variance of the indices are also derived. A significant advantage of the proposed framework is that only the parametrized Markov chain generator matrix is required as an input, i.e., closed-form expressions for the reliability and performability indices as a function of the model parameters are not needed. Several case studies illustrate the accuracy of the proposed method in approximating distributions of reliability and performability indices. Additionally, analysis of a large model demonstrates lower execution times compared to Monte Carlo simulations.