Revision history [back]

This is not a solution but only a very awkward attempt which illustrates what is missing (coefficient_list is home-brewed).

egf = 1 + sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
tay = integral(taylor(egf, x, 0, 6),x)
r = coefficient_list(tay, x, true)
red = lambda i,o: o/i if i>0 else o
P = [red(ind, obj) for ind, obj in enumerate(r[1::])]  
for narayana_poly in P:
    print coefficient_list(narayana_poly, t)  

[1]
[0, 1]
[0, 1, 1]
[0, 1, 3, 1]
[0, 1, 6, 6, 1]
[0, 1, 10, 20, 10, 1]
[0, 1, 15, 50, 50, 15, 1]

I suspect Sage will not get many friends at the OEIS if it cannot cope with such elementary maneuvers with generating functions.