You can use polynomial curve fitting (eqivalent to using Lagrange Interpolating Polynomials) to solve it:
Recall that there is precisely one polynomial of degree n that passes through n+1 distinct points. It's simple in concept, but the formula is somewhat ugly, where: Suppose you are given three points, say
\left( x_{1}, y_{1}\right) ,\mbox{ } \left( x_{2}, y_{2}\right) ,\mbox{ and } \left( x_{3}, y_{3}\right)
then the polynomial that interpolates (that is fits) those points is of degree 3-1=2, and will consist of three terms. It is
P(x) = y_{1}\frac{\left( x-x_{2}\right) \left( x-x_{3}\right) }{\left( x_{1}-x_{2}\right) \left( x_{1}-x_{3}\right)} + y_{2}\frac{\left( x-x_{1}\right) \left( x-x_{3}\right) }{\left( x_{2}-x_{1}\right) \left( x_{2}-x_{3}\right)} + y_{3}\frac{\left( x-x_{1}\right) \left( x-x_{2}\right) }{\left( x_{3}-x_{1}\right) \left( x_{3}-x_{2}\right)}
since we wanted it to have the value y_{1} at x_{1}, the first term does this and it is zero at x_{2}\mbox{ and } x_{3} (this is the numerator) and the denominator ensures that at x_{1} the numerator terms are corrected for so that it is still y_{1} there. The other terms likewise.
To sum the k-th powers of the first n integers first note that the sum is a polynomial in n of degree k+1 (you have proved this for k=1,...,5 already), which can thus be completely specified by k+2 points, say the first k+2 points. Then find the polynomial in terms of n for (each k) by polynomial curve fitting as such:
\sum_{i=1}^{n} i^{k}=\sum_{i=1}^{k+2}\left(\sum_{m=1}^{i} m^{k} \prod_{j=1}^{k+2} \frac{n-j}{i-j} \right)
where the terms in the product corresponding to j=i are omitted.
The \sum_{m=1}^{i} m^{k} for i=1,..., k+2 terms are the y_{i}\mbox{'s} , the "the first k+2 points" mentioned above.
I hope this was clear enough for you to follow. There's also a way using the Bernoulli numbers, I wasn't fond of it though.
-Ben