$n$th derivative of $e^{1/x}$
I am trying to find the $n$'th derivative of $f(x)=e^{1/x}$. When looking at the first few derivatives I noticed a pattern and eventually found the following formula
$$\frac{\mathrm d^n}{\mathrm dx^n}f(x)=(-1)^n e^{1/x} \cdot \sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 n+k}$$
I tested it for the first $20$ derivatives and it got them all. Mathematica says that it is some hypergeometric distribution but I don't want to use that. Now I am trying to verify it by induction but my algebra is not good enough to do the induction step.
Here is what I tried for the induction (incomplete, maybe incorrect)
$\begin{align*} \frac{\mathrm d^{n+1}}{\mathrm dx^{n+1}}f(x)&=\frac{\mathrm d}{\mathrm dx}(-1)^n e^{1/x} \cdot \sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 n+k}\\ &=(-1)^n e^{1/x} \cdot \left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} (-2n+k) x^{-2 n+k-1}\right)-e^{1/x} \cdot \sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 (n+1)+k}\\ &=(-1)^n e^{1/x} \cdot \sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k}((-2n+k) x^{-2 n+k-1}-x^{-2 (n+1)+k)})\\ &=(-1)^{n+1} e^{1/x} \cdot \sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k}(2n x-k x+1) x^{-2 (n+1)+k} \end{align*}$
I don't know how to get on from here.
How's this?
$$\left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} \left(-2n+k\right) x^{-2 n+k-1}\right) - \left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 \left(n+1\right)+k}\right) =$$
$$= \left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} \left(-2n+k\right) x^{-2\left(n+1\right)+k+1}\right) - \left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 \left(n+1\right)+k}\right) =$$
$$= \left(\sum _{k'=1}^{n} \left(k'-1\right)! \binom{n}{k'-1} \binom{n-1}{k'-1} \left(-2n+k'-1\right) x^{-2\left(n+1\right)+k'}\right) - \left(\sum _{k=0}^{n-1} k! \binom{n}{k} \binom{n-1}{k} x^{-2 \left(n+1\right)+k}\right) =$$
$$= \left(\sum _{k'=0}^{n} \left(k'-1\right)! \binom{n}{k'-1} \binom{n-1}{k'-1} \left(-2n+k'-1\right) x^{-2\left(n+1\right)+k'}\right) - \left(\sum _{k=0}^{n} k! \binom{n}{k} \binom{n-1}{k} x^{-2 \left(n+1\right)+k}\right) =$$
$$= \sum _{k=0}^{n} \left(\left(k-1\right)! \binom{n}{k-1} \binom{n-1}{k-1} \left(-2n+k-1\right) - k! \binom{n}{k} \binom{n-1}{k}\right) x^{-2 \left(n+1\right)+k}$$
Then $$\left(k-1\right)! \binom{n}{k-1} \binom{n-1}{k-1} \left(-2n+k-1\right) - k! \binom{n}{k} \binom{n-1}{k} =$$
$$= \frac{\left(k-1\right)!n!\left(n-1\right)!\left(-2n+k-1\right)}{\left(n-k+1\right)!\left(k-1\right)!\left(n-k\right)!\left(k-1\right)!} - \frac{k!n!\left(n-1\right)!}{\left(n-k\right)!k!k!\left(n-k-1\right)!} =$$
$$= \frac{n!\left(n-1\right)!\left(-2n+k-1\right)k}{\left(n-k+1\right)!\left(n-k\right)!k!} - \frac{n!\left(n-1\right)!\left(n-k\right)\left(n-k+1\right)}{\left(n-k\right)!k!\left(n-k+1\right)!} =$$
$$= \frac{n!\left(n-1\right)!}{\left(n-k+1\right)!\left(n-k\right)!k!} \left(\left(-2n+k-1\right)k - \left(n-k\right)\left(n-k+1\right)\right) =$$
$$= \frac{-n\left(n+1\right)n!\left(n-1\right)!}{\left(n-k+1\right)!\left(n-k\right)!k!} =$$
$$= -\frac{\left(n+1\right)!}{\left(n-k+1\right)!} \binom{n}{k} =$$
$$= -k! \binom{n+1}{k} \binom{n}{k}$$
Peter's and Mike's answers have clearly settled this question; I'll just explain the OP's mention of "Mathematica says that it is some hypergeometric distribution". More specifically, one wonders how Mathematica might have arrived at the Kummer confluent hypergeometric function ${}_1 F_1\left({{a}\atop{b}}\mid x\right)$.
We start with the transformed Maclaurin series:
$$\exp\left(\frac1{x}\right)=1+\sum_{k=1}^\infty \frac1{x^k k!}$$
and then differentiate the series termwise, using the formula
$$\frac{\mathrm d^n}{\mathrm dx^n}x^{-k}=\frac{(-1)^n (n+k-1)!}{(k-1)! x^{k+n}}$$
to yield
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n}{x^n}\sum_{k=1}^\infty \frac{(n+k-1)!}{k!(k-1)! x^k}$$
which can be reindexed like so:
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n}{x^n}\sum_{k=0}^\infty \frac{(n+k)!}{(k+1)!k! x^{k+1}}$$
and rewritten in terms of Pochhammer symbols to yield
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n}{x^n}\sum_{k=0}^\infty \frac{n!(n+1)_k}{(2)_k k! x^{k+1}}$$
or, expressed hypergeometrically,
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n n!}{x^{n+1}}{}_1 F_1\left({{n+1}\atop{2}}\mid \frac1{x}\right)$$
It's not immediately obvious that this one can be expressed in terms of $\exp(1/x)$ multiplied by a finite sum, but we can then apply the Kummer transformation to the hypergeometric function like so:
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n n!}{x^{n+1}}\exp\left(\frac1{x}\right){}_1 F_1\left({{1-n}\atop{2}}\mid -\frac1{x}\right)$$
and find that the hypergeometric function has a nonpositive numerator parameter $1-n$ for $n > 0$, which means it can be expressed as a terminating series (since the Pochhammer symbol $(1-n)_k$ is $0$ for $k \geq n$):
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n n!}{x^{n+1}}\exp\left(\frac1{x}\right)\sum_{k=0}^{n-1}\frac{(-1)^k (1-n)_k}{k!(k+1)!x^k}$$
Rewriting the Pochhammer symbol in terms of factorials, we get
$$\frac{\mathrm d^n}{\mathrm dx^n}\exp\left(\frac1{x}\right)=\frac{(-1)^n n!}{x^{n+1}}\exp\left(\frac1{x}\right)\sum_{k=0}^{n-1}\frac{(n-1)!}{(n-k-1)!k!(k+1)!x^k}$$
which is equivalent to the formula you obtained.
My point here, I guess, is that hypergeometric functions are nothing to be uncomfortable with; if it helps, just consider them as shorthand for certain (in)finite sums.
Here's a completely elementary argument that I like better than my previous answer that uses the Faà di Bruno formula. All that is required are the Maclaurin series for $e^x$ and the definition of the Lah numbers $L(n,k)$ as the coefficients arising when converting rising factorial powers to falling factorial powers via $$x^{\overline{n}} = \sum_{k=1}^n L(n,k) x^{\underline{k}}.$$
First, substitute $1/x$ for $x$ in the Maclaurin series for $e^x$ and differentiate $n$ times. This yields $$\frac{d^n}{dx^n} e^{1/x} = (-1)^n x^{-n} \sum_{j=0}^{\infty} \frac{j^{\overline{n}} x^{-j}}{j!}.$$
Then convert rising powers to falling powers to obtain
$$\frac{d^n}{dx^n} e^{1/x} = (-1)^n x^{-n} \sum_{j=0}^{\infty} \sum_{k=1}^n L(n,k) \frac{j^{\underline{k}} x^{-j}}{j!}.$$
Working with the summation in $j$, we have $$\sum_{j=0}^{\infty} \frac{j^{\underline{k}} x^{-j}}{j!} = \sum_{j=k}^{\infty} \frac{j^{\underline{k}} x^{-j}}{j!} = \sum_{j=k}^{\infty} \frac{j! x^{-j}}{(j-k)! j!} = \sum_{i=0}^{\infty} \frac{x^{-i-k}}{i!} = x^{-k} e^{1/x},$$ where the last step uses, once again, the Maclaurin series for $e^x$.
Thus we have $$\frac{d^n}{dx^n} e^{1/x} = (-1)^n e^{1/x} \sum_{k=1}^n L(n,k) x^{-n-k},$$ which is just your formula with an index change, as $$L(n,k) = \binom{n}{k}\binom{n-1}{k-1} (n-k)!.$$
(From another perspective, this is computing the rising factorial moments of a Poisson $(1/x)$ distribution. With this in mind, converting to falling factorial powers seems reasonable, as the falling factorial moments of a Poisson distribution have a known simple expression. See also this related answer to a previous math.SE question.)
You can get the same result using the Faà di Bruno formula and Bell polynomials.
The Faà di Bruno formula says that $$\frac{d^n}{dx^n} f(g(x)) = \sum_{k=0}^n f^{(k)}(g(x)) B_{n,k}\left(g'(x),g''(x),\dots,g^{(n-k+1)}(x)\right),$$
where $B_{n,k}(x_1,x_2,\dots,x_{n-k+1})$ is a partial Bell polynomial.
Since $f(x) = e^x$, and $g(x) = 1/x$, we have $$\frac{d^n}{dx^n} e^{1/x} = \sum_{k=0}^n e^{1/x} B_{n,k}\left(-x^{-2},2x^{-3},\dots,(-1)^{n-k+1}(n-k+1)!x^{-n+k-2}\right).$$
It is known that $B_{n,k}(1!,2!,\dots,(n-k+1)!) = \binom{n}{k}\binom{n-1}{k-1} (n-k)!$, a Lah number (see Wikipedia on Bell polynomials or Comtet's Advanced Combinatorics, p. 135). Thus (keeping track of signs and powers of $x$), $$\frac{d^n}{dx^n} e^{1/x} = (-1)^n e^{1/x} \sum_{k=0}^n \binom{n}{k}\binom{n-1}{k-1} (n-k)! x^{-n-k},$$ which is just your formula with an index change (and where $\binom{n-1}{-1} = 0$).