Integral and integro-differential equations have many applications in various fields of science and engineering such as biological models, industrial mathematics, control theory of financial mathematics, economics, fluid dynamics, heat and mass transfer, queuing theory, electrostatics, electromagnetic, electrodynamics, elasticity, biomechanics, oscillation theory, and so forth. It is well known that it is extremely difficult to analytically solve nonlinear integro-differential equations. Indeed, few of these equations can be solved explicitly. So it is required to devise an efficient approximation scheme for solving these equations. So far, several numerical methods are developed. The solution of the first order integro-differential equations has been obtained by the numerical integration methods such as Euler-Chebyshev and Runge-Kutta methods. Moreover, a differential transform method for solving integro-differential equations was introduced by Darania and Ebadian. Shidfar et al. applied the homotopy analysis method for solving the nonlinear Volterra and Fredholm integro-differential equations. As a concrete example, we can express the mathematical model of cell-to-cell spread of HIV-1 in tissue cultures considered by Mittler et al. Babolian suggested an effective direct method to determine the numerical solution of the specific nonlinear Volterra-Fredholm integro-differential equations. Their approach was based on triangular functions. In recent years, the meshless methods have gained more attention not only by mathematicians but also in the engineering community. An efficient hybrid method is developed to approximate the solution of the high-order nonlinear Volterra-Fredholm integro-differential equations. The proposed approximation scheme is based on a hybrid of block-pulse functions and the well known Lagrange interpolating polynomials. The hybrid of block pulse functions and Lagrange interpolating polynomials was originally introduced by Marzban and Hoseini. They implemented this kind of hybrid functions for solving Volterra’s population model. Because block- pulse functions and Lagrange interpolating polynomials are both complete, therefore, the proposed hybrid functions constitutes a complete set in the Hilbert space L 2 [0 ,t f ). The nice properties of hybrid functions together with the associated operational matrices of integration, product and derivative with the Kronecker property of these functions are then used to reduce the solution of the nonlinear Volterra- Fredholm integro-differential equations to the solution of a nonlinear system of algebraic equations whose solution is much more easier than the original one. Various problems arising in science and engineering can be modeled by nonlinear Volterra-Fredholm integral equations, the main focus of this study is to present an effective numerical method for solving them. The mentioned operational matrices are sparse, hence making the method computationally attractive without sacrificing the accuracy of the solution. Moreover, It is shown that a few number of hybrid functions are needed to achieve high accuracy and a satisfactory convergence. The validity and applicability of the proposed method are demonstrated through illustrative examples. The method is general, easy to implement and yields the desired accuracy in a few terms of hybrid basis functions. The simulation results demonstrate the effectiveness of the suggested numerical scheme.