If the matrix function f(At) posses the properties of f(At)=gf(tkA, then the recurrence formula fi−1=gfi,i=N,N−1,⋯,1,f(tA)=f0, can be established. Here, fN=f(AN)=∑j=0majANj,AN=tkNA. This provides an algorithm for computing the matrix function f(At). By specifying the calculation accuracy p, a method is presented to determine m and N in a way that minimizes the time of the above algorithm, thus providing a fast algorithm for f(At). It is important to note that m only depends on the calculation accuracy p and is independent of the matrix A and t. Therefore, f(AN) has a fixed calculation format that is easily computed. On the other hand, N depends not only on A, but also on t. This provides a means to select t such that N is equal to 0, a property of significance. In summary, the algorithm proposed in this article enables users to establish a desired level of accuracy and then utilize it to select the appropriate values for m and N to minimize computation time. This approach ensures that both accuracy and efficiency are addressed concurrently. We develop a general algorithm, then apply it to the exponential, trigonometric, and logarithmic matrix functions, and compare the performance with that of the internal system functions of Mathematica and Pade approximation. In the last section, an example is provided to illustrate the rapid computation of numerical solutions for linear differential equations.