<p>This study introduces an efficient numerical scheme for solving the fractional-order FitzHugh-Nagumo (FN) equation under given initial and boundary conditions. Unlike traditional methods that often pair high-order spatial discretization with low-order temporal schemes, our approach ensures a balanced, high-order treatment of both time and space. The method employs fractional Lagrange basis functions, constructed on Chebyshev-Gauss-Lobatto collocation points, to interpolate the solution. This formulation transforms the original fractional differential problem into a nonlinear algebraic system via an operator matrix, which is then solved using the Levenberg-Marquardt algorithm. A rigorous convergence analysis is provided. The method’s superior accuracy and efficiency are demonstrated through six illustrative examples, with favorable comparisons made against several recent numerical techniques. The numerical simulations also yield biological insights: we examine how the neuronal excitation threshold influences the action potential, and we model epileptic neuronal behavior by incorporating a physiologically relevant input function.</p>