Note: This answer is a supplement to the beautiful answer of @Winther. It's in no way an alternative, but an additional examination. OPs question which is R.P. Stanleys problem $47c$ in Enumerative Combinatorics Ed.02 and his elaboration give rise to some questions and I was curious to find some information.
We know from @Winther that due to the commutator property of the operator $(x+D)$ a special case of the BCH Theorem can be applied leading to
\begin{align*}
e^{(x+D)t}=e^{xt}e^{Dt}e^{\frac{t^2}{2}}\tag{1}
\end{align*}
Interesting aspects:
(1) OPs approach calculating $(x+D)^n$ iteratively for small $n$ leads to matrices with a relation to binomial coefficients in some skew diagonals . What is the connection with (1) and how could we derive a solution based upon OPs approach?
(2) Since R.P. Stanleys problems $47a$ and $47b$ are closed related to the current one and they show identities with Stirling numbers $s(n,k)$ of the first kind and Stirling numbers $S(n,k)$ of the second kind we should expect, that $47c$ is also related to Stirling numbers.
(3) The first row of OPs matrix for $n=5$ OPs indicates a relation to Hermite polynomial. How do they interact?
Aspect (1) OPs matrices and Bessel numbers
We introduce for convenience only the following shorthand notation
\begin{align*}
(x+D^\ast)^n:=\sum_{k=0}^n\binom{n}{k}x^{n-k}D^k
\end{align*}
Note: This is not an operator in the common sense, since we don't have e.g. $(x+D^\ast)^{n+1}=(x+D^\ast)(x+D^\ast)^{n}$.
Setting $n=4$ and $n=5$ we obtain a representation of $(x+D)^n$ as sum of binomial expressions of the form $(x+D^\ast)^k$
\begin{align*}
(x+D)^4&=x^4+4x^3D+6x^2D^2+4xD^3+D^4\\
&\qquad+6(x^2+2xD+D^2)+3\\
&=\color{blue}{\mathbf{1}}(x+D^\ast)^4+\color{blue}{\mathbf{6}}(x+D^\ast)^2+\color{blue}{\mathbf{3}}
\\
\\
(x+D)^5&=x^5+5x^4D+10x^3D^2+10x^2D^3+5xD^4+D^5\\
&\qquad+10(x^3+3x^2D+3xD^2+D^3)+15(x+D)\\
&=\color{blue}{\mathbf{1}}(x+D^\ast)^5+\color{blue}{\mathbf{10}}(x+D^\ast)^3+\color{blue}{\mathbf{15}}(x+D^\ast)
\end{align*}
This representation shows the binomial coefficients and the multiples we can see in the diagonals of OPs matrices. We focus on the coefficients of $(x+D^\ast)^k$ and find
\begin{align*}
\begin{array}{rrrrrrrrrr}
&k&0&1&2&3&4&5&6&7\\\hline
n&&\phantom{105}&\phantom{105}&\phantom{105}&\phantom{105}
&\phantom{105}&\phantom{105}&\phantom{105}&\phantom{105}\\
0&&1\\
1&&&1\\
2&&&1&1\\
3&&&&3&1\\
4&&&&\color{blue}{\mathbf{3}}&\color{blue}{\mathbf{6}}&\color{blue}{\mathbf{1}}\\
5&&&&&\color{blue}{\mathbf{15}}&\color{blue}{\mathbf{10}}&\color{blue}{\mathbf{1}}\\
6&&&&&15&45&15&1\\
7&&&&&&105&105&21&1\\
\end{array}
\end{align*}
The blue entries show the coefficients of the calculation above. Now OEIS identifies these entries as sequence A100861: Triangle of Bessel numbers read by rows.
The following information is based upon The Bessel numbers and matrices by Sheng Liang Yang and Zhan Ke Qiao and
Bessel numbers are related to Stirling Numbers of the Second kind, since:
Bessel numbers $B(n,k)$ are the number of partitions of an $n$-set into $k$ non empty subsets, each of size at most $2$.
Stirling numbers of the second kind $S(n,k)$ are the number of partitions of an $n$-set into $k$ non empty subsets.
The connection can also nicely be seen via the vertical exponential generating functions for $B(n,k)$ and $S(n,k)$
\begin{align*}
\sum_{n=k}^{\infty}B(n,k)\frac{t^n}{n!}&=\frac{1}{k!}\left(t+\frac{t^2}{2}\right)^k\tag{2}\\
\\
\sum_{n=k}^{\infty}S(n,k)\frac{t^n}{n!}&=\frac{1}{k!}\left(e^t-1\right)^k\\
&=\frac{1}{k!}\left(t+\frac{t^2}{2!}+\frac{t^3}{3!}+\ldots\right)^k
\end{align*}
The number of summands in the right hand side of (2) indicates the restriction of the size of the $k$ subsets by at most $2$.
We now obtain a representation of $(x+D)^n$ with Bessel numbers:
\begin{align*}
(x+D)^n&=B(n,n)(x+D^\ast)^n+B(n,n-1)(x+D^\ast)^{n-1}\\
&\qquad+\ldots
+B(n,n-\left\lfloor\frac{n}{2}\right\rfloor)(x+D^\ast)^{n-2\left\lfloor\frac{n}{2}\right\rfloor}\\
&=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}B(n,n-k)(x+D^\ast)^{n-2k}\\%\tag{3}\\
&=\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}B(n,n-k)\sum_{l=0}^{n-2l}\binom{n-2k}{l}x^lD^{n-2k-l}\\
\end{align*}
The coefficients of $(x+D)^n=\sum_{i,j}a_{n,i,j}x^iD^j$ are in terms of Bessel numbers
\begin{align*}
a_{n,i,j}&=[x^iD^j](x+D)^n\\
&=[D^j]\sum_{k=0}^{\lfloor\frac{n}{2}\rfloor}B(n,n-k)\sum_{l=0}^{n-2l}\binom{n-2k}{i}D^{n-2k-i}\\
&=\begin{cases}
{\displaystyle B(n,\frac{n+i+j}{2})\binom{i+j}{i}}&\qquad\qquad i+j\equiv n(2)
\\
0&\qquad\qquad \text{otherwise}
\end{cases}
\end{align*}
In the last line we use $n-2k-i=j$, so $k=\frac{n-i-j}{2}$ if $i+j\equiv n(2)$
Connection with $e^{xt}e^{Dt}e^{\frac{t^2}{2}}$:
Using (1) we obtain following representation of $B(n,k)$.
\begin{align*}
e^{xt}e^{Dt}e^{\frac{t^2}{2}}&=\sum_{i\geq 0}\frac{x^i}{i!}t^i \sum_{j\geq 0}\frac{D^j}{j!}t^je^{\frac{t^2}{2}}\\
&=\sum_{i\geq 0}\sum_{j\geq 0}\binom{i+j}{i}x^iD^j\frac{t^{i+j}}{(i+j)!}\sum_{k\geq 0}\frac{1}{2^k}\frac{t^{2k}}{k!}\\
&=\sum_{n\geq 0}\sum_{{i,j,k\geq 0}\atop{i+j+2k=n}}\binom{i+j}{i}\frac{1}{(i+j)!}\frac{1}{k!}\frac{1}{2^k}x^iD^jt^{i+j+2k}\\
&=\sum_{n\geq 0}\sum_{{0\leq i,j\leq n}\atop{i+j\equiv n(2)}}\binom{i+j}{i}\frac{n!}{(i+j)!\left(\frac{n-i-j}{2}\right)!}\cdot\frac{1}{2^{\frac{n-i-j}{2}}}x^iD^j\frac{t^{n}}{n!}\\
\end{align*}
We conclude
\begin{align*}
B\left(n,\frac{n+i+j}{2}\right)=
\begin{cases}
{\displaystyle
\frac{n!}{(i+j)!\left(\frac{n-i-j}{2}\right)!}\cdot\frac{1}{2^{\frac{n-i-j}{2}}}
}&\qquad\qquad i+j\equiv n(2)\\
0&\qquad\qquad \text{otherwise}
\end{cases}
\end{align*}
Of course, we could also derive a closed formula for $B(n,k)$ with other means, e.g. the representation of (2).
$$$$
Aspect (2) Connection with Stirling numbers
A well known relationship between Stirling numbers of the first kind and of the second kind is
\begin{align*}
\sum_{j=0}^ns(n,j)S(j,k)=\delta_{n,k}\qquad n\geq 0, \quad 0\leq k \leq n
\end{align*}
with $\delta_{n,k}$ the Kronecker delta symbol. The authors use exponential Riordan arrays in their paper to show in (19) a similar identity with Bessel numbers.
\begin{align*}
B(n,k)=\sum_{j=k}^n2^{j-k}s(n,j)S(j,k)\qquad n\geq 0, \quad 0\leq k \leq n
\end{align*}
$$$$
Aspect (3) Connection with Hermite Polynomials
We can also find information regarding Hermite Polynomials in the paper. The sequence of row sums of the table with the Bessel numbers above can according to Theorem 3.2 be represented by Hermite polynomials,
\begin{align*}
H_n(x)=\sum_{k=0}^{\left\lfloor\frac{n}{2}\right\rfloor}\frac{(-1)^kn!}{k!(n-2k)!}(2x)^{n-2k}
\end{align*}
which can be defined by the generating function
\begin{align*}
\sum_{n=0}^{\infty}H_n(x)\frac{t^n}{n!}=e^{2xt-t^2}
\end{align*}
The shape of this generating function makes a connection with the generating functions of the Bessel numbers plausible. In fact, the Hermite functions play an important role when studying differentiation operators. See e.g. Generalized polynomials, operational identities and their applications by G. Dattoli. He presents with formula (18) the following generalisation of $(x+D)^n$ in terms of Hermite polynomials:
\begin{align*}
(\alpha D+\beta x)^n=\sum_{k=0}^n\binom{n}{k}\alpha^kH_{n-k}\left(\beta x,\frac{1}{2}\alpha\beta\right)D^k
\end{align*}
