Yesterday I struggled with the Narayana polynomials which can be computed with Sage (and the help of kcrisman) as

```
P = lambda n: hypergeometric([-n,-n+1],[2], 1/x).simplify_hypergeometric()
[expand(x^k*P(k)) for k in (0..7)]
```

Today I want to generate the coefficients of these polynomials by an exponential generating function. With Maple I have no problems.

```
egf := 1 + int((sqrt(t)*exp((1+t)*x)*BesselI(1,2*sqrt(t)*x))/x,x);
s := n -> n!*coeff(series(egf,x,n+2),x,n);
seq(print(seq(coeff(s(n),t,j),j=0..n)),n=0..6);
```

This gives the requested triangle:

```
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
```

With Sage I tried:

```
from sage.symbolic.integration.integral import indefinite_integral
t = var('t')
h = lambda x, t: sqrt(t)*exp((1+t)*x)*bessel_I(1,2*sqrt(t)*x)
egf = 1 + indefinite_integral(h(x,t),x)
taylor(egf, x, 0, 6)
```

This gives me a TypeError. Perhaps this is the wrong Ansatz.

Any help is appreciated.