series seems to make errors . luckily, taylor works correctly.
start with L_F=(exp(-s)-1+s)/(s^2/2) (laplace transform of the "equilibrium" uniform density)
and switch to Pollaczek laplace transform L_F /(1- epsilon * L_F);
You can check that the moments expansions t, t1 differ
var('s');n=2 L_F=(exp(-s)-1+s)/(s^2/2) #satisfies L_F(s=0)=1 L_L=L_F/(1-L_F/3) t = L_L.series(s,2n+2) t1= taylor(L_L,s,0,2n+2) print t1