Phase type (PH) distributions are widely used in modeling and simulation due to their generality and analytical properties. In such settings, it is often necessary to construct a PH distribution that aligns with real-world data by matching a set of prescribed moments. Existing approaches provide either exact closed-form solutions or iterative procedures that may yield exact or approximate results. However, these methods are limited to matching a small number of moments using PH distributions with a small number of phases, or are restricted to narrow subclasses within the PH family. We address the problem of approximately fitting a larger set of given moments using potentially large PH distributions. We introduce an optimization methodology that relies on a re-parametrization of the Markovian representation, formulated in a space that enables unconstrained optimization of the moment-matching objective. This reformulation allows us to scale to significantly larger PH distributions and capture higher moments. Results on a large and diverse set of moment targets show that the proposed method is, in the vast majority of cases, capable of fitting as many as 20 moments to PH distributions with as many as 100 phases, with small relative errors on the order of under 0.5% from each target. We further demonstrate an application of the optimization framework where we search for a PH distribution that conforms not only to a given set of moments but also to a given shape. Finally, we illustrate the practical utility of this approach through a queueing application, presenting a case study that examines the influence of the i^{th} moment of the inter-arrival and service time distributions on the steady-state probabilities of the GI/GI/1 queue length.