In this paper, we present an operational matrix method to obtain numerical solutions of 𝜓-Caputo fractional ordinary and partial differential equations. For this purpose, a fractional version of the Taylor theorem is presented in the framework of 𝜓-fractional calculus. The method converts the underlying ordinary or partial differential equations to systems of algebraic equations. The method is accompanied by numerical examples to verify the applicability and effectiveness of the proposed method. Further, estimates of upper bounds of error for the approximations have been derived.