I found the following formula in my previous question. This differs from my previous question in that I want an alternative proof of the below recursive formula for calculating $\displaystyle\sum_{k=1}^nk^p$.
Suppose I had a function recursively defined as
$$f(x,p)=a_px+p\int_0^xf(t,p-1)dt$$
$$a_p=1-p\int_0^1f(t,p-1)dt$$
For $p\in\mathbb N$. For $p=0$, we trivially get $f(x,0)=x$, which shall be our initial condition.
It can then be noticed that
$$a_1=1-\int_0^1tdt=\frac12$$
$$f(x,1)=\frac12x+\int_0^xtdt=\frac12x+\frac12x^2$$
$$a_2=1-2\int_0^1\frac12t+\frac12t^2dt=\frac16$$
$$f(x,2)=\frac16x+\int_0^x\frac12t+\frac12t^2dt=\frac16x+\frac14x^2+\frac16x^3$$
And the general pattern is $f(x,p)=\sum_{k=1}^xk^p$ whenever $x\in\mathbb N\quad(?)$
How do I prove that whenever $x\in\mathbb N$
$$f(x,p)=\sum_{k=1}^xk^p$$
without applying the methods mentioned in the link above?
Define $$B_p(x):=p\left(f(x,p-1)-x^{p-1}\right)+a_p$$ then $$B'_p(x)=p\left(f'(x,p-1)-(p-1)x^{p-2}\right)=p\left [\left(a_{p-1}x+(p-1)\int_0^xf(t,p-2)\textrm{d}t\right)^{'}-(p-1)x^{p-2} \right]=p\left(a_{p-1}+(p-1)f(x,p-2)-(p-1)x^{p-2}\right)=pB_{p-1}(x)$$ and $$\int_0^1B_p(t)\textrm{d}t=\int_0^1\left[p\left(f(t,p-1)-t^{p-1}\right)+a_p\right]\textrm{d}t=p\int_0^1f(t,p-1)\textrm{d}t-p\int_0^1t^{p-1}\textrm{d}t+a_p=1-a_p-p\frac{1}{p}+a_p=0$$ With the starting value $B_0(w)=1$ this uniquely determines the sequence of Bernoulli polynomials see Koenigsberger, K. (2003) Analysis 1: Springer.
It follows that $$B_p(0)=p\left(f(0,p-1)-0^{p-1}\right)+a_p=a_p=B_p$$ and $$\frac{B_{p+1}(x+1)-B_{p+1}}{p+1}=\frac{(p+1)(f(x+1,p)-(x+1)^p)+a_{p+1}-a_{p+1}}{p+1}=f(x+1,p)-(x+1)^p=\sum_{k=1}^xk^p.$$