With the gradual transitioning to Network function virtualization (NFV), Internet service providers face the problem of where to introduce NFV in order to make the most benefit of that; here, we measure the benefit by the amount of traffic that can be served in an NFV-enabled network. This problem is non-trivial as it is composed of two challenging subproblems: 1) placement of nodes to support virtual network functions (referred to as VNF-nodes); 2) allocation of the VNF-nodes' resources to network flows. This problem has been studied for the one-dimensional setting, where all network flows require one network function, which requires a unit of resource to process a unit of flow. In this work, we consider the multi-dimensional setting, which introduces new challenges in addition to those of the one-dimensional setting (e.g., NP-hardness and non-submodularity). To address these difficulties, we propose a novel two-level relaxation method that allows us to draw a connection to the sequence submodular theory and utilize the property of sequence submodularity along with the primal-dual technique to design two approximation algorithms. We prove that the proposed algorithms have a non-trivial approximation ratio and perform trace-driven simulations to show the effectiveness of the proposed algorithms.