The Taylor operational matrix of integration and the Taylor product operational matrix are introduced. These two operational matrices are applied to approximation of solutions of Fredholm and Volterra integral equations. The proposed method reduces solution of integral equations to the successive solution of a set of linear algebraic equations in matrix form. Owing to the simplicity of the operational matrix of integration, and the product operational matrix of the Taylor series, the algorithms derived possess considerable computational advantages over the orthogonal-polynomial approximation, provided that both input and output are analytic functions oft.