restart:
Coefficient in a polynomial
with(gfun):
P := (1+t)^10000 * (1+t + t^3)^60000;
factor(diff(P, t)/P);
diffeq := {diff(y(t), t) - %*y(t), y(0)=1};
rec := diffeqtorec(diffeq, y(t), u(n));
pr := rectoproc(rec, u(n)):
pr(1000);
# coeff(expand(P), t, 1000); # ARGH!!
coeff(series(P, t, 1001),t, 1000); # OK, but slow
Guessing
L := [seq(add(binomial(n,k)^3, k=0..n), n=0..12)];
f := add(L[i+1]*t^i, i=0..nops(L)-1);
add(add(a[i,j]*t^i*`if`(j > 0, t^j*diff(f, t$j), f), i=0..2), j=0..2);
series(%, t, 13);
eqs := [seq(coeff(%, t, i), i=0..12)];
solve(%);
subs(%, add(add(a[i,j]*t^(i+j)*D^j, i=0..2), j=0..2));
collect(%, D, factor);
Impemented in gfun
listtodiffeq(L, y(t))[1];
rec := SumTools[Hypergeometric][ZeilbergerRecurrence](binomial(n,k)^3, n, k, u, 0..n);
rectodiffeq({rec, u(0)=1, u(1)=2}, u(n), y(t));