The derivative of Γ(z) is given by Ψ(z)Γ(z), where Ψ(z) is the digamma function. Higher derivatives are given in terms of Ψn(z), the polygamma function.

For z = 1, Ψn(1) = (−1)n+1n!ζ(n + 1), where ζ(z) is the Riemann zeta function. For large n, ζ(n) is asymptotically 1.

From this, you should be able to compute the coefficients explicitly.

When I look at the coefficients numerically, they do not appear to converge towards 0 that rapidly. They do become smaller eventually, but the progression is uneven:

1, -.577, -.656, .420, .167, .0422, -.00962, -.00722, -.00117, .000215, .000128, ...

Thank you very much.

Actually the first 26 coefficients of a similar function 1/Gamma(z) are under formula 6.1.34 of Abramowitz Stegun "Pocketbook of Mathematical Functions"

The cancellations that lead to the fast convergence of the coefficients are crucial and are not obvious from the first couple of coefficients as you rightfully write.For  cancellations the high digit precision is very important.

So I have calculated the coefficients to high precision (>100 digits) up to order n=100:

For faster convervegence I have used

exp(EULER*z)/Gamma(1-z) = exp(sum(-zeta(n)/n*z^n))  (can be deduced from 6.1.33 of the above mentioned book.

I have found my formula from these numerical calculations with a coefficient A = 2/PI or at least quite close to it.

But I had hoped someone had done a saddlepoint calculation of the coefficients for high order (hopefully that is possible) and so a good formula was already known somewhere in the literature. 

Both Γ(1-z) and  1/Γ(1-z) are analytic in the neighbourhood of z=0. The sought-after asymptotic series for 1/Γ(1-z) for z → 0 is thus obtained from the Taylor series expansion of Γ(1-z) around z=0 followed by the calculation of the reciprocal of this series. For the calculation of the former series, make use of the expression in Eq. (6.1.1) of Abramowitz & Stegun; because of the exponential function exp(-t) in the integrand of this expression, the coefficients of any finite-order Taylor series expansion of Γ(1-z) around z=0 are bounded.

Incidentally, if you have Mathematica at hand, calculation of the sought-after asymptotic series is very straightforward.*  With γ = 0.577215... denoting the Euler constant, one has:

1//Γ(1-z) ~ 1 - γ z + (6 γ2 - π2) z2/12 + ... 

______

* Just use Series[1/Gamma[1-z], {z,0,n}], where n is a given finite integer. For the n-th coefficient in this series, with n a given integer, use SeriesCoefficient[1/Gamma[1-z],{z,0,n}].

Thank you very much Behnam.

But my problem is not to find a way to calculate the first couple of Taylor-coefficients 1/Gamma(1-z) = sum(a_n*z^n) analytically or numerically. My problem is to find the asymptotic form of a_n for large n.

I have calculated the first 100 coefficients numerically (and of course the first few coefficients a_0 = 1, a_1 = -gamma, a_2 = (6*gamma^2 - PI^2)/12 analytically) and as you write this is quite straightforward.

But how does a_n behave in leading order for big n? As I wrote numerically

a_n ~ 1/(n!)^A*cos(B*n + C)*n^D*E (1 + o(1)).

seems to be feasible with A = 2/PI or at least close.

Large orders of Taylor coefficients can of course often be calculated by saddle point evaluations. In our case because there are infinitely many saddle points it looks difficult though. So I still hope someone has solved this before.

I think I finally found the answer in Theorem 3.1. of http://arxiv.org/abs/1407.5983 by Lazhar Fekih-Ahmed. So the formula is more involved than I thought and have written above and contains the Lambert W function. But as usual it is found by saddle point evaluation as I had thought.