With the large-scale integration of renewable energy and power electronic devices into power grids, harmonic issues in distribution networks have become increasingly complex due to multi-source interactions, posing significant challenges to traditional modeling methods. This paper proposes a harmonic source modeling method based on the Bernstein-Kolmogorov-Arnold Network (Bernstein-KAN), aiming to address the high training costs and computational complexity of existing machine learning models when handling multiple harmonic features. By integrating the Kolmogorov-Arnold theorem, Bernstein-KAN replaces traditional neural network weight parameters with Bernstein polynomials, enhancing model flexibility and interpretability. Additionally, an Adaptive Factor Model (AFM) is introduced for dimensionality reduction to further optimize computational efficiency. Experimental validation using measured harmonic data compared Bernstein-KAN with models such as MLP, Spline-KAN, and Cheby-KAN. The results demonstrate that Bernstein-KAN achieves an average Mean Absolute Percentage Error (MAPE) of 4.714% for 3rd, 5th, 7th, and 11th harmonic predictions, with a training time of only 0.312 s, significantly outperforming other methods. Compared to the AFM-MOFGPR model, Bernstein-KAN reduces training time by over 99% in scenarios with multiple harmonic features while maintaining prediction accuracy. The study concludes that Bernstein-KAN substantially lowers computational costs without compromising precision, offering an efficient solution for harmonic source modeling in complex power grid environments.