I would like to find the series expansion of
(E^(x^k/k!) Gamma[k, x])/Gamma[k]
for $k$ being a positive integer, up to the order of $x^{2k+1}$. Mathematica provides a function Series, which can deal with this when $k$ is replaced by any fixed integers. But it cannot handle the case $k$ remains a symbol.
Is it actually possible to get this for general $k$ in Mathematica?
If I correctly understand your question, this can be done in such a way:
Expand[Normal[Series[Exp[y], {y, 0, 2}]] /. y -> x^k/k!*
Normal[Series[Gamma[k, x]/Gamma[k], {x, 0, 2},
Assumptions -> k \[Element] Integers && k > 0]]]
$$\frac{x^{4 k}}{2 k^2 (k!)^2 \Gamma (k)^2}+\frac{x^{2 k}}{2 (k!)^2}-\frac{x^{2 k}}{k k! \Gamma (k)}-\frac{x^{3 k}}{k (k!)^2 \Gamma (k)}+\frac{x^{2 k+1}}{(k+1) k! \Gamma (k)}-\frac{x^{2 k+2}}{2 (k+2) k! \Gamma (k)}+\frac{x^{3 k+1}}{(k+1) (k!)^2 \Gamma (k)}-\frac{x^{3 k+2}}{2 (k+2) (k!)^2 \Gamma (k)}-\frac{x^{4 k+1}}{k (k+1) (k!)^2 \Gamma (k)^2}+\frac{x^{4 k+2}}{2 (k+1)^2 (k!)^2 \Gamma (k)^2}+\frac{x^{4 k+2}}{2 k (k+2) (k!)^2 \Gamma (k)^2}-\frac{x^{4 k+3}}{2 (k+1) (k+2) (k!)^2 \Gamma (k)^2}+\frac{x^{4 k+4}}{8 (k+2)^2 (k!)^2 \Gamma (k)^2}+\frac{x^k}{k!}+1 $$
Increasing the order of the power series expansion (two in the above), one obtains a huge expression.
Here is a slightly different method for doing it:
ode = First[Head[DifferentialRootReduce[GammaRegularized[k, x], x]]];
ser = Normal[Series[Exp[z], {z, 0, 3}]] /. z -> x^k/k!;
ser2 = FullSimplify[AsymptoticDSolveValue[ode[y, x], y[x], {x, 0, 3}]];
Assuming[k > 0, Expand[FullSimplify[FunctionExpand[ser ser2]]]]
1 + x^(1 + k)/((1 + k) Gamma[k]) - x^(2 + k)/(2 (2 + k) Gamma[k]) +
x^(3 + k)/(6 (3 + k) Gamma[k]) - x^(4 k)/(6 Gamma[1 + k]^4) +
(k x^(1 + 4 k))/(6 (1 + k) Gamma[1 + k]^4) -
(k x^(2 + 4 k))/(12 (2 + k) Gamma[1 + k]^4) +
(k x^(3 + 4 k))/(36 (3 + k) Gamma[1 + k]^4) - x^(3 k)/(3 Gamma[1 + k]^3) +
(k x^(1 + 3 k))/(2 (1 + k) Gamma[1 + k]^3) -
(k x^(2 + 3 k))/(4 (2 + k) Gamma[1 + k]^3) +
(k x^(3 + 3 k))/(12 (3 + k) Gamma[1 + k]^3) - x^(2 k)/(2 Gamma[1 + k]^2) +
(k x^(1 + 2 k))/((1 + k) Gamma[1 + k]^2) -
(k x^(2 + 2 k))/(2 (2 + k) Gamma[1 + k]^2) +
(k x^(3 + 2 k))/(6 (3 + k) Gamma[1 + k]^2)
E^(x^k/k!) GammaRegularized[k, x]. $\endgroup$