In solving a mathematical problem numerically, we frequently need to operate on a vector υ by an operator that can be expressed asf(A), whereA is anN ×N matrix [e.g., exp(A), sin(A), A−-]. Except for very simple matrices, it is impractical to construct the matrixf (A) explicitly. Usually an approximation to it is used. This paper develops an algorithm based upon a polynomial approximation tof (A). First the problem is reduced to a problem of approximatingf (z) by a polynomial in z, where z belongs to a domainD in the complex plane that includes all the eigenvalues ofA. This approximation problem is treated by interpolatingf (z) in a certain set of points that is known to have some maximal properties. The approximation thus achieved is “almost best.” Implementing the algorithm to some practical problems is described. Since a solution to a linear systemAx=b isx=A−1b, an iterative solution algorithm can be based upon a polynomial approximation tof (A)=A−1. We give special attention to this important problem.