Markov reward models have interesting modeling applications, particularly those addressing fault-tolerant hardware/software systems. In this paper, we consider a Markov reward model with a reward structure including only reward rates associated with states, in which both positive and negative reward rates are present and null reward rates are allowed, and develop a numerical method to compute the distribution function of the cumulative reward till exit of a subset of transient states of the model. The method combines a model transformation step with the solution of the transformed model using a randomization construction with two randomization rates. The method introduces a truncation error, but that error is strictly bounded from above by a user-specified error control parameter. Further, the method is numerically stable and takes advantage of the sparsity of the infinitesimal generator of the transformed model. Using a Markov reward model of a fault-tolerant hardware/software system, we illustrate the application of the method and analyze its computational cost. Also, we compare the computational cost of the method with that of the (only) previously available method for the problem. Our numerical experiments seem to indicate that the new method can be efficient and that for medium-size and large models can be substantially faster than the previously available method.