This article presents a procedure for computing the renewal function, the renewal density, the integral of the renewal function, and the variance function of phase-type renewal processes. The procedure hinges on the computation of the state probability vector of a continuous-time Markov chain. This is accomplished by using a randomization approach that is simple, efficient, and numerically stable and does not require numerical integration. I discuss approximating arbitrary interrenewal distributions by phase-type distributions so that the procedure can be used to approximate renewal and related functions.