Let σ0,σ1,⋯,σn be a set of n+ 1 distinct real numbers (i.e., σi≠σj, for i≠j) and F0,F1,⋯ ,Fn, be given real s × r matrices, we know that there exists a unique s × r matrix polynomial Pn(λ) of degree n such that Pn(σi) = Fi, for i = 0,1,⋯ ,n, Pn is the matrix interpolation polynomial for the set {(σi,Fi),i = 0,1,⋯ ,n}. The matrix polynomial Pn(λ) can be computed by using the Lagrange formula or the barycentric method. This paper presents a new method for computing matrix interpolation polynomials. We will reformulate the Lagrange matrix interpolation polynomial problem and give a new algorithm for giving the solution of this problem, the Recursive Matrix Polynomial Interpolation Algorithm (RMPIA) in full and simplified versions, and some properties of this algorithm will be studied. Cost and storage of this algorithm with the classical formulas will be studied and some examples will also be given.