The Kuramoto model (KM) of coupled phase oscillators on graphs provides the most influential framework for studying collective dynamics and synchronization. It exhibits a rich repertoire of dynamical regimes. Since the work of Strogatz and Mirollo [J. Stat. Phys., 63 (1991), pp. 613-635], the mean field equation derived in the limit as the number of oscillators in the KM goes to infinity has been the key to understanding a number of interesting effects, including the onset of synchronization and chimera states. In this work, we study the mathematical basis of the mean field equation as an approximation of the discrete KM. Specifically, we extend the Neunzert's method of rigorous justification of the mean field equation (cf. [H. Neunzert, Fluid Dyn. Trans., 9 (1978), pp. 229-254]) to cover interacting dynamical systems on graphs. We then apply it to the KM on convergent graph sequences with non-Lipschitz limit. This family of graphs includes many graphs that are of interest in applications, e.g., nearest-neighbor and small-world graphs. The approaches for justifying the mean field limit for the KM proposed previously in [C. Lancellotti, Transp. Theory Statist. Phys., 34 (2005), pp. 523-535; H. Chiba and G. S. Medvedev, arXiv:1612.06493, 2016] do not cover the non-Lipschitz case.