An incorrect method to sum the first $n$ squares which nevertheless works
Solution 1:
Let $f_k(n)=\sum_{i=1}^n i^k$. We all know that $f_k$ is actually a polynomial of degree $k+1$. Also $f_k$ can be characterised by the two conditions: $$f_k(x)-f_k(x-1)=x^k$$ and $$f_k(0)=0.$$ Differentiating the first condition gives $$f_k'(x)-f_k'(x-1)=k x^{k-1}.$$ Therefore the polynomial $(1/k)f_k'$ satisfies the first of the two conditions that $f_{k-1}$ does. But it may not satisfy the second. But then $(1/k)(f_k'(x)-f_k'(0))$ does. So $$f_{k-1}(x)=\frac{f_k'(x)-f_k'(0)}k.$$
The mysterious numbers $f_k'(0)$ are related to the Bernoulli numbers, and when $k\ge3$ is odd they obligingly vanish...
Solution 2:
Consider Bernoulli's formula which gives
$$ \sum_1^n i^m = \frac{1}{m+1} \sum_0^m {m+1 \choose k} B_k n^{m+1-k} $$
where $B_k$ is the k'th Bernoulli number. Take the derivative with respect to $n$ on the right hand side, you get
$$ \frac{d}{dn} \sum_1^n i^m = \frac{1}{m+1} \sum_0^m \frac{(m+1)\cdot m!}{(m+1 - k)\cdot (m -1 + 1 - k)! k!}\cdot (m+1-k) B_k n^{m-1 + 1-k} $$
$$= \sum_0^m { m \choose k} B_k n^{m-k} = m \sum_1^n i^{m-1} + B_m$$
So the difference between formally taking the derivative w.r.t. $i$ as you did, and formally taking the derivative w.r.t. $n$, is exactly the Bernoulli constant for the power $m$. Now Bernoulli's constant is equal to 0 for all odd $m$ that is greater than 1. For those cases your observation is also true: the two "derivatives" are equal.
Now, some random guessing, since I am not an analytic number theorist. The deep connection to Bernoulli numbers probabaly has something to do with analytic number theory and harmonic analysis that I do not know myself. That the formulae work for odd powers is possibly something related to the method of descent, where analytical results in even dimensions can be arrived from results in one higher dimension via some sort of trace.
Solution 3:
This is quite obvious if expressed most naturally in terms of the Bernoulli polynomials, viz.
$\quad\quad\quad\quad\rm\ \ \ \sum n^k\ = \:\frac{1}{k+1} (B_{k+1}(n+1) - B_{k+1}(0))\ $
and $\rm\quad\quad\ B_k^{\:'}(x)\ =\ k\ B_{k-1}(x)\ $
thus $\ \rm\ (\sum n^k)'\ =\ k\ \sum n^{k-1} + B_{\:k} $
According to Knuth's exposition these identities go all the way back to Jacobi - with special cases known to Faulhaber. Note that the identities in Robin Chapman's post are simply special cases of these well-known identities for Bernoulli polynomials.
The Bernoulli polynomials are a special-case of polynomials amenable to study by the powerful techniques of the Umbral Calculus. This is the best way to understand the genesis of their many interesting properties. For a nice introduction see Steven Roman's book.