The present study presents an efficient numerical method for solving the fractional differential equations governing electrohydrodynamic flows. This method is based on an effective combination of Bernoulli polynomial bases and hyperbolic sine functions. In this regard, a collocation method is developed that employs combined Bernoulli–hyperbolic bases for approximation. First, the highest-order fractional derivative appearing in the equation is approximated through an expansion based on fractional-order Bernoulli polynomials and the hyperbolic sine function. Then, by applying the Riemann–Liouville fractional integral operator and imposing the boundary conditions of the problem, approximations for the solution function as well as its lower-order fractional derivatives are obtained. Finally, by forming the residual function and collocating at the nodes corresponding to the roots of Legendre polynomials, a system of nonlinear algebraic equations is derived, which is solved using Newton’s algorithm to obtain the numerical solution of the problem. Also, convergence analysis of the proposed method is presented and evaluation results show that the proposed method achieves accuracy with a small number of collocation points.