[20-08] I will add something that sheds light on why the latter works out nicely. In what follows we provide a proof of Faulhaber's Formula using the operation of binomial convolution on sequences of real numbers (in fact, we needn't worry about the whole real numbers, but merely the rational numbers). This operation arises naturally when multiplying (formal) exponential generating series of the form $$\sum_{n\geqslant 0}\gamma_n \frac{z^n}{n!}$$ and collecting powers of $z$. The problem of finding $1^n+2^n+\cdots+m^n$ motivates us to look into each of the sequences $1,j,j^2,j^3,\ldots$ that have the EGF $e^{jz}$. Thus the sum has EGF $$\sum_{j=1}^m e^{jz}=e^z\frac{e^{mz}-1}{e^{z}-1}=\frac{e^{mz}-1}z \frac{(-z)}{e^{-z}-1}$$
It is easy to obtain explicit generating sequences for $\dfrac{e^{mz}-1}z$; indeed $$\frac{e^{mz}-1}z=\sum_{n\geqslant 0}\frac{m^{n+1}}{n+1}\frac{z^n}{n!}$$
In particular when $m=1$ we get the Harmonic numbers generate $$H(z)=\frac{e^{z}-1}z$$
The Bernoulli numbers arise naturally as the inverses of the Harmonic numbers, in the sense that
$$B(z)=\frac{z}{e^z-1}=\sum_{n\geqslant 0}B_n\frac{z^n}{n!}$$
The first equation then gives Faulhaber's Formula immediately: we have $$\sum_{j=1}^m e^{jz}=e^z\frac{e^{mz}-1}{e^{z}-1}=B(-z)\frac{e^{mz}-1}z $$
Equating coefficients and using the convolution gives $$\sum_{j=1}^m j^n =\sum_{k=0}^n\binom nk (-1)^k B_k\frac{m^{n-k+1}}{n-k+1}$$
This is the shortest proof I know of. We work things formally, without worrying about convergence. We use the fact that a sequence with non-zero first term $a_0$ admits a unique inverse (as is the case of the Harmonic numbers). This is proven by observing the convolution gives a recurrence for the terms of this inverse.
The above is very similar to what will follow. The first lemma gives $$H(z)e^{jz}=\frac{e^{(j+1)z}-1}z-\frac{e^{jz}-1}z$$ and the second gives $$e^z\frac{e^{mz}-1}z=\frac{e^{(m+1)z}-1}z-\frac{e^{z}-1}z$$ Telescoping gives $$H(z)\sum\limits_{j = 1}^m {{e^{jz}}} = \frac{{{e^{(m + 1)z}} - {e^z}}}{z} = {e^z}\frac{{{e^{mz}} - 1}}{z}$$
Thus $$\sum\limits_{j = 1}^m {{e^{jz}}} = \frac{{{e^{(m + 1)z}} - {e^z}}}{z} = {e^z}B\left( z \right)\frac{{{e^{mz}} - 1}}{z}$$ and using $B(-z)=e^zB(z)$ $-$ this $1\star B=\mu B$ which was quite ugly to prove without EGFs $-$ finishes things to give Faulhaber's Formula $$\sum\limits_{j = 1}^m {{e^{jz}}} = B\left( { - z} \right)\frac{{{e^{mz}} - 1}}{z}$$
DEF Let ${\bf x} =(x_0,x_1,x_2,\ldots,x_n,\ldots)$, ${\bf y}=(y_0,y_1,y_2,\ldots,y_n,\ldots)$ be sequences. We define a new sequence, their binomial convolution, by the formula $$\left({\bf x}\star {\bf y}\right)_n:=\sum_{k=0}^n \binom{n}{k}x_ky_{n-k}$$
OBS The binomial convolution is associative and commutative.
DEF We define the special sequences $\mu=(1,-1,1,-1,1,-1,1,\ldots)$, ${\bf 1 }=(1,1,1,1,1,\ldots)$, ${\rm id}=(1,0,0,0,0,0,\ldots)$. Note then that $$\mu\star {\bf 1} =\rm{id}$$ $${\rm id}\star {\bf x}={\bf x}$$
We're ready to prove the
THM (Inversion formula) Let ${\bf x} =(x_0,x_1,x_2,\ldots,x_n,\ldots)$, ${\bf y}=(y_0,y_1,y_2,\ldots,y_n,\ldots)$ be sequences. Then ${\bf x}\star {\bf 1}={\bf y}$ if, and only if ${\bf x}={\bf y}\star \mu$.
P By the above $$\begin{align}{\bf x}\star {\bf 1}&={\bf y}\\({\bf x}\star {\bf 1})\star \mu&={\bf y}\star \mu\\{\bf x}\star ({\bf 1}\star\mu)&={\bf y}\star \mu\\{\bf x}\star {\rm id }&={\bf y}\star \mu\\{\bf x}&={\bf y}\star \mu\end{align}$$
and the other direction is analogous.
[7-28] I have decided to define the Bernoulli numbers as the inverse of the Harmonic numbers, since this puts in evidence their importance when proving Faulhaber's formula. I will keep looking for a proof of the closed form formula.
DEF Denote by $\bf H$ the Harmonic numbers $\left(1,\frac 1 2,\frac 13,\frac 14 \ldots\right)$. We define the Bernoulli numbers by $${\bf H}\star {\bf B}=\rm id$$
This is well-defined since inverses are unique. A few values are $$\left(1,-\dfrac 1 2,\dfrac 1 6,0,-\dfrac 1 {30},0,\dfrac1{42},0,-\dfrac{1}{30},\ldots\right)$$
OBS The definition says that $\sum_{k=0}^n\binom{n+1}kB_k=[n=0]$ and ${\bf B}\star {\bf 1}={\bf B}+[n=1]$
PROP For $n\geqslant 1$ we have $B_{2n+1}=0$.
P (Credits to Rob) From ${\bf B}\star {\bf 1}={\bf B}+[n=1]$ we obtain by inversion that ${\bf B} ={\bf B}\star \mu+[n=1]\star \mu$. But $[n=1]\star \mu=-(-1)^nn$; so ${\bf B}\star \mu=B+\mu{\bf N}$. Since ${\bf x}\star \mu {\bf y}=\mu(\mu{\bf x}\star {\bf y})$ we obtain $\mu{\bf B}\star {\bf 1}=\mu{\bf B}+{\bf N}$. Now we write things explicitly $$\sum_{k=0}^n\binom{n+1}{k}B_k=[n=0]\;\;,\;\;\sum_{k=0}^{m-1}\binom mk(-1)^kB_k=m$$
Thus by $m\mapsto n+1$ $$\sum_{k=0}^n\binom{n+1}{k}B_k=[n=0]\;\;,\;\;\sum_{k=0}^{n}\binom {n+1}k(-1)^kB_k=n+1$$
Thus $$n+1-[n=0]=\sum_{k=0}^n\left((-1)^k-1\right)\binom{n+1}kB_k$$ and $n\mapsto 2n+1$ gives $$2n+2=\sum_{k=0}^{2n+1}\left((-1)^k-1\right)\binom{2n+2}kB_k$$ but the $k=1$ term is $-2\binom{2n+2}1\left(-\frac1 2\right)=2n+2$ so, since even terms vanish we get $$0=-2\sum_{k=1}^{n}\binom{2n+2}{2k+1}B_{2k+1}$$ and induction does the rest. $\blacktriangle$
COR The formula ${\bf B}\star{\bf 1}=\mu{\bf B}$ holds.
First, we have:
LEMMA1 Fix $j\geqslant 0$. Then, $$\left({\bf{H}} \star {j^k}\right)_n = \frac{{{{\left( {j + 1} \right)}^{n + 1}} - {j^{n + 1}}}}{{n + 1}}$$
P This follows from direct calculation and use of the binomial theorem.
LEMMA2 Fix $m$. Then $$\left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right)_n = \frac{{{{\left( {m + 1} \right)}^{n + 1}} - 1}}{{n + 1}}$$
P This should also be fairly obvious by integrating the binomial expansion and finding the appropriate constant of integration.
COR By the two previous lemmas: $$\left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{{\left( {{\bf{H}} \star {j^k}} \right)}_n}} $$
Now we can prove
THM (Faulhaber's Formula) $$\left( {\mu {\bf{B}} \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{j^n}} $$
P We have that $$\sum\limits_{j = 1}^m {{{\left( {{\bf{H}} \star {j^k}} \right)}_n}} = \left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right)$$
By the distributive property $$\left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = {\bf{H}} \star \sum\limits_{j = 1}^m {{j^n}} $$ and after multiplication by $\bf B$, this is $$\eqalign{
& {\bf{B}} \star \left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{j^n}} \cr
& \left( {\mu {\bf{B}} \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{j^n}} \cr} $$
as desired. $\blacktriangle$
Best Answer
A simple-minded approach is to solve the two variable recurrence relation iteratively, that is, knowing $a_{n,0}$ find $a_{n,1}$, then $a_{n,2}$, etc.
We must have
$$\begin{eqnarray*} a_{n,1} &=& a_{n-1,1}+a_{n-1,0} \\ &=& a_{n-1,1}+1, \qquad a_{1,1}=1. \end{eqnarray*}$$ This is a one variable recurrence relation of the form $$b_n = b_{n-1}+1, \qquad b_1 = 1.$$ This can be solved by the usual methods. We find $a_{n,1} = n.$
Next we have $$\begin{eqnarray*} a_{n,2} &=& a_{n-1,2}+a_{n-1,1} \\ &=& a_{n-1,2}+n-1, \qquad a_{2,2}=1. \end{eqnarray*}$$ This is another simple recurrence relation. We find $a_{n,2} = n(n-1)/2.$
At the next step, $$\begin{eqnarray*} a_{n,3} &=& a_{n-1,3}+a_{n-1,2} \\ &=& a_{n-1,3} + \frac{1}{2}(n-1)(n-2), \qquad a_{3,3}=1. \end{eqnarray*}$$ This implies $a_{n,3} = n(n-1)(n-2)/6.$
This process can be repeated to build up $a_{n,k}$ for any $k$. At some point a pattern will be noticed and the principle of induction can be applied.