Method 1: Higher-order Pearson systems
The Pearson system is, by convention, taken to be the family of solutions $p(x)$ to the differential equation:
$$ \frac{d p (x)}{dx} \; = \; -\frac{(a+x) }{c_0 + c_1 x + c_2 x^2} \; p(x) $$
where the four Pearson parameters $(a, c_0, c_1, c_2)$ can be expressed in terms of the first four moments of the population.
Instead of basing the Pearson system upon the quadratic $c_0 + c_1 x + c_2 x^2$, one can instead consider using higher order polynomials as the foundation stone. So, for example, one can consider a Pearson-style system based upon a cubic polynomial. This will be the family of solutions $p(x)$ to the differential equation:
$$ \frac{d p(x)}{dx} \; = \; -\frac{(a+x) }{c_0 + c_1 x + c_2 x^2 + c_3 x^3} \; p (x) $$
which yields the solution:
I solved this for fun some time back (having the same thought train as the OP): the derivation and solution is given in Chapter 5 of our book; if interested, a free download is available here:
http://www.mathstatica.com/book/bookcontents.html
Note that whereas the second-order (quadratic) Pearson family can be expressed in terms of the first 4 moments, the third-order (cubic) Pearson-style family requires the first 6 moments.
Method 2: Gram-Charlier expansions
Gram-Charlier expansions are also discussed in the same Chapter 5 (see section 5.4) ... and also allow one to construct a fitted density, based on arbitrarily large $k^{th}$ moments. As the OP suggests, the Gram-Charlier expansion expresses the fitted pdf as a function of a series of derivatives of the standard Normal pdf, known as Hermite polynomials. The Gram-Charlier coefficients are solved as a function of the population moments ... and the bigger the expansion, the more moments required. You may also wish to look at related Edgeworth expansions.
Population moments or sample moments??
For the Pearson-style system: If the moments of the population are known, then using higher moments should unambiguously yield a better fit. If, however, the observed data is a random sample drawn from the population, there is a trade-off: a higher order polynomial implies that higher order moments are required, and the estimates of the latter may be unreliable (have high variance), unless the sample size is 'large'. In other words, given sample data, fitting using higher moments can become 'unstable' and produce inferior results. The same is true for Gram-Charlier expansions: adding an extra term can actually yield a worse fit, so some care is required.