I have been trying to find the sequence $a_0,a_1,a_2,\dots$ that solves the following nested recurrence relation: $$ \begin{aligned} (1) &&c_0 &=a_0^m\\ (2) &&c_k &=\frac{1}{a_0k}\sum_{\ell=1}^k((m+1)\ell-k)a_\ell c_{k-\ell}\\ (3) &&0 &=(1+2m)c_k-m(c_{k-1}+\sum_{\ell=0}^{k-1}a_\ell c_{k-\ell-1}), \end{aligned} $$ where $m=2,4,6,\dots$ is some positive even integer. Solving this system directly proved intractable so I tried generating the first several $a_k$'s and then made a guess. I came up with $$ a_k=\frac{m^{k-1}}{(1+2m)^k k!}\sum_{\ell=0}^k\binom{k}{\ell}(1+\tfrac{\ell+1}{m})_{k-1}a_0^{\ell+1}, $$ where $(s)_n=\Gamma(s+n)/\Gamma(s)$ is the Pochhammer symbol (rising factorial). Indeed, this guess holds for $k=0,1,\dots,15$. My code for verifying this guess gets really slow for $k>20$ and so this is as far as I could compare my guess with the exact values for $a_k$.
My question is how to prove my guess is correct (or incorrect) for all $k\in\Bbb N$. It would seem I must first substitute my guess into (2) and solve for $c_k$. Then I could see if the solution for $c_k$ solves (3). However, solving for $c_k$ also seems intractable which has me stuck. I figured there must be an easier way to verify if this guess holds or not. Can I verify without explicitly solving for $c_k$?
Edit:
Here is my Mathematica code:
a[k_]:= Sum[m^(k-1)/((1+2 m)^k k!) Binomial[k,l] Pochhammer[1+(l+1)/m,k-1]Subscript[a,0]^(l+1),{l,0,k}];
c[0]:=Subscript[a,0]^m;
c[k_]:=1/(Subscript[a,0] k) Sum[((m+1) l-k) a[l] c[k-l],{l,1,k}];
Eqn[k_]:=(1+2 m) c[k]-m (c[k-1]+Sum[a[l] c[k-l-1],{l,0,k-1}]);
Table[FullSimplify[Eqn[k],m>=2 && Element[m/2,Integers]],{k,1,15,1}]
Context:
The $a_k$'s are the coefficients of a perturbation series expansion for the real $x$ that solves $$ \frac{\mathrm d}{\mathrm dx}(1+x+x^2+\dots+x^m)=0. $$
Update:
Using Faa di Bruno's formula I was able to write a "closed-form" for $c_k$. We have, $$ \begin{aligned} c_0 &=a_0^m\\ c_k &=\frac{1}{k!}\sum_{\ell=1}^k(m)^{(\ell)}a_0^{m-l}B_{k,\ell}(1!\, a_1,\dots,(k-\ell+1)!\,a_{k-\ell+1}), \end{aligned} $$ where $(s)^{(n)}=\Gamma(s+1)/\Gamma(s-n+1)$ is the falling factorial and $B_{n,k}(x_1,\dots,x_{n-k+1})$ is the partial Bell-polynomial. We can implement this in Mathematica as
a[l_]:= Sum[(Subscript[a,0] m^(l-1))/((1+2 m)^l l!) Binomial[l,j] Pochhammer[1+(j+1)/m,l-1] Subscript[a,0]^j,{j,0,l}];
c[0]:=Subscript[a,0]^m;
c[k_]:=1/k! Sum[FactorialPower[m,l] Subscript[a,0]^(m-l) BellY[k,l,Table[j! a[j],{j,1,k-l+1}]],{l,1,k}];