In this paper, we present an efficient method for computing the action of matrix $\phi$-functions. Our approach is based on a scaling and recursing procedure and incorporates a shifting technique as a preprocessing step to enhance efficiency. We conduct a forward error analysis to determine the optimal scaling parameter and polynomial degree for achieving the desired accuracy. Numerical comparisons with existing algorithms demonstrate that the proposed algorithm performs well in terms of both accuracy and efficiency.