Abstract: Meshfree collocation methods are favored owing to their inherent simplicity in computer implementation and high efficiency in numerical computation. However, numerical results indicate that the order of kernel functions in conjunction with support sizes significantly affects the accuracy of meshfree collocation methods, as underscores the critical need for choosing an optimized kernel function in order to improve the accuracy of meshfree collocation computation. In meshfree collocation analysis, the kernel functions usually are chosen by trial and error, and a theoretical way to optimize the kernel functions is still missed. In this study, through the accuracy analysis for meshfree collocation methods, a quantitative method is proposed for optimally selecting kernel functions, where a direct link is established between the kernel function parameters and the accuracy of numerical solutions. In order to further improve the efficiency, the second order gradients of meshfree shape functions are constructed using the implicit first order gradients, which avoids the costly computation of second-order meshfree gradients and then improves the computational efficiency. Subsequently, a local truncation error analysis leads to an error expression for meshfree collocation methods. It turns out that the error coefficient can be extracted from the error expression and serve as an index that effectively measures the relationship between the kernel functions and the accuracy of meshfree collocation methods. Consequently, an optimal kernel function is determined via minimizing the error coefficients across different kernel orders and support sizes for meshfree collocation computation. It is noteworthy that the proposed error coefficient for meshfree collocation methods does not require solving any system of equations, and is quite general and independent on specific problems. The effectiveness of the proposed quantitative method regarding the determination of optimal kernel functions for meshfree collocation methods is well demonstrated through the excellent agreement between the theoretically predicted optimal kernel functions and the corresponding accuracy superiority of numerical results.