Behaviors of many dynamic systems with uncertainty can be modelled effectively by systems of fuzzy differential equations. In this paper, we develop new numerical iterative method for solving systems of fuzzy initial value problems based on the reproducing kernel theory under the assumption of Hukuhara differentiability. The exact and approximate solutions are given with series form in terms of their parametric form, where two smooth reproducing kernel functions are used throughout the evolution of the algorithm to obtain the required nodal values. Furthermore, error estimation is proved in order to capture the behavior of fuzzy solutions. Applicability, potentiality, and efficiency of the proposed algorithm for the fuzzy solutions of different fuzzy systems are investigated using computer tables and graphical representation. © 2017 NSP Natural Sciences Publishing Cor.