Network Function Virtualization (NFV) on Software-Defined Networks (SDN) can effectively optimize the allocation of Virtual Network Functions (VNFs) and the routing of network flows simultaneously. Nevertheless, most previous studies on NFV focus on unicast service chains and thereby are not scalable to support a large number of destinations in multicast. On the other hand, the allocation of VNFs has not been supported in the current SDN multicast routing algorithms. In this paper, therefore, we make the first attempt to tackle a new challenging problem for finding a service forest with multiple service trees, where each tree contains multiple VNFs required by each destination. Specifically, we formulate a new optimization, named Service Overlay Forest (SOF), to minimize the total cost of all allocated VNFs and all multicast trees in the forest. We design a new $3rho_{ST}$-approximation algorithm to solve the problem, where $rho_{ST}$ denotes the best approximation ratio of the Steiner Tree problem, and the distributed implementation of the algorithm is also presented. Simulation results on real networks for data centers manifest that the proposed algorithm outperforms the existing ones by over 25%. Moreover, the implementation of an experimental SDN with HP OpenFlow switches indicates that SOF can significantly improve the QoE of the Youtube service.