<p>This paper proposes a new and efficient numerical method for solving time-fractional partial differential equations. This approach, which we call the fractional Bernoulli-Picard iteration method, takes advantage of the strengths of both Bernoulli functions and the well-established Picard iteration method. Picard’s method is a popular iterative algorithm for solving initial value problems. However, it is challenging to calculate integrals involving complex and non-linear functions. Our proposed method addresses this limitation by introducing fractional-order Bernoulli functions as an approximation tool in the integral. A significant advantage of fractional Bernoulli functions is the presence of a tunable parameter that can be optimized for high accuracy. Notably, the fractional integrals of Bernoulli functions can be easily calculated at each iteration step and significantly simplify the computational burden. To strengthen the theoretical foundation of our method, we carefully analyze its convergence properties. We further verify its efficiency and accuracy through numerical simulations and demonstrate its potential as a powerful tool in solving practical problems including fractional-order convection-diffusion, wave-like, and inhomogeneous Burgers equations.</p>