This paper presents a new numerical method for solving fractional optimal control problems (FOCPs). The fractional derivative in the dynamic system is described in the Caputo sense. The method is based upon Bernoulli polynomials. The operational matrices of fractional Riemann–Liouville integration and multiplication for Bernoulli polynomials are derived. The error upper bound for the operational matrix of the fractional integration is also given. The properties of Bernoulli polynomials are utilized to reduce the given optimization problems to the system of algebraic equations. By using Newton’s iterative method, this system is solved and the solution of FOCPs are achieved. Illustrative examples are included to demonstrate the validity and applicability of the technique.