This issue can be treated using 2D convolution.
- First of all, instead of the initial hexagon, we are going to work in an equivalent but distorted one (see Fig. 1) in order that the points of interest (the endpoints of paths with $n$ "steps") have simple coordinates ; in this way, one has just to build an evolving matrix representation (see program below). The "distorded" hexagon is represented on Fig. 1 for the case $n=2$ :

Fig. 1. Case $n=2$. Left figure : endpoints are materialized by the number of "paths" reaching them. These points have coordinates of the form $$\left(i/2^{n+1},j/2^{n+1}\right)$$ where $(i,j)$ are integers indexing the entries of the attached matrix. Right figure : a simple affine transform (a transvection) gives back a "classical hexagon" ; further (simple) affine transforms are needed to transform this hexagon into your one $A=\left(1,0\right), B=\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right), $ etc. But a straightforward solution is obtained by using barycentric coordinates (see in particular fig. 3 and the associated Edit below).
- One can have an idea of the process by looking at Fig. 2 representing this convolution operation by "unfolding" the first steps in a third dimension (zeros are ommited):

Fig. 2 : The presence of convolution can be understood through this 3D representation of the issue, reminding neural networks. We have taken here a normalization such that all the $(x,y)$ coordinates are in $[0,1] \times [0,1]$.
Based on this idea, one can write the following very compact program (written in Matlab but translatable into any other scientific language) :
S=[1,1,0;
1,0,1;
0,1,1];
M=[1];
ke=3;
for k=1:ke;
M=kron(M,[1,0;0,0]);
M=conv2(M,S);
end;
Comments on this program :
The $6$ "ones" in matrix $S$ account for the oblique directions in which the information is "spread", replacing the six directions $(\cos(k\pi/3),\sin(k\pi/3)$. (Letter $S$ reminds the strong analogy with the so-called "structuring element" in the domain of "mathematical morphology").
The first line inside the for-loop accounts for the progressive refinement of the working mesh (approximate doubling of the number of lines and columns at each new step) : more precisely, between each pair of lines (resp. columns) of the matrix at rank $n$, a line (resp. column) of zeros is introduced. This operation is realized by using a Kronecker product.
The second line operates a 2D convolution.
Optionaly, one can add just before closing the "for loop", the following instruction
p=2^(k+1)-1;M=M(1:p,1:p)
"cropping" matrix $M$ in order to cancel (harmless...) trailing zeros on the right and bottom margins.
- the end value $ke$ of $k$ can be anything... sensible for your computer.
The program above gives instantly
the successive arrays :
$$\begin{array}{|ccc|}\hline
1&1&0\\
1&0&1\\
0&1&1\\
\hline\end{array}$$
$$\begin{array}{|ccccccc|}\hline
1&1&1&1&0&0&0\\
1&0&2&0&1&0&0\\
1&2&1&1&2&1&0\\
1&0&1&0&1&0&1\\
0&1&2&1&1&2&1\\
0&0&1&0&2&0&1\\
0&0&0&1&1&1&1\\
\hline\end{array}$$
$$\begin{array}{|cccccccccccccc|}\hline 1&1&1&1&1&1&1&1&0&0&0&0&0&0&0\\
1&0&2&0&2&0&2&0&1&0&0&0&0&0&0\\
1&2&1&1&3&3&1&1&2&1&0&0&0&0&0\\
1&0&1&0&2&0&2&0&1&0&1&0&0&0&0\\
1&2&3&2&1&3&3&1&2&3&2&1&0&0&0\\
1&0&3&0&3&0&2&0&3&0&3&0&1&0&0\\
1&2&1&2&3&2&1&1&2&3&2&1&2&1&0\\
1&0&1&0&1&0&1&0&1&0&1&0&1&0&1\\
0&1&2&1&2&3&2&1&1&2&3&2&1&2&1\\
0&0&1&0&3&0&3&0&2&0&3&0&3&0&1\\
0&0&0&1&2&3&2&1&3&3&1&2&3&2&1\\
0&0&0&0&1&0&1&0&2&0&2&0&1&0&1\\
0&0&0&0&0&1&2&1&1&3&3&1&1&2&1\\
0&0&0&0&0&0&1&0&2&0&2&0&2&0&1\\
0&0&0&0&0&0&0&1&1&1&1&1&1&1&1\\
\hline \end{array}$$
etc.
Edit (2023/08/28) : Have a look at Fig. 3 and its notations.

Fig. 3.
It is obtained (see program below) by the use of barycentrical coordinates in triangle $GHI$.
The idea here is that indexes $(i,j)$ of matrix $M$ can be used for the generation of barycentrical coordinates, under the condition to be conveniently "resized". Let us explain how.
Let us begin by index $j$ (the column index of matrix $M$) running from $1$ to $p=2^n-1$. $j/a$ plays the role of "degree of attraction" of vertex $G$ for the current point with normalization coefficient $a:=(3/2)(2^n)$. Why ? Because barycentric coordinates of points inside hexagon $ABCDEF$ with respect to triangle $GHI$ have $2/3$ as their extreme value ;
Dealing with index $i$ (the row index), as it is decreasing, it is a "degree of repulsion" with respect to vertex $H$. In order to get a "degree of attraction", we have to consider $(2^n-i)/a$
the last degree of attraction (with repect to vertex $I$), its formula is taken in such a way that the sum of the 3 degrees of attraction (i.e., the barycentrical coordinates) of the current point is equal to $1$.
Fig. 3 is plainly generated (set apart the names of vertices and the drawings of the hexagon and the triangle) by placing the following lines at the bottom of the first program :
n=ke+1;d=2^n;p=d-1;a=(3/2)*d;
s=sqrt(3);G=[3/2;-s/2];H=[0;s];I=[-3/2;-s/2];
for i=1:p
for j=1:p
V=(j*G+(d-i)*H+(a+i-d-j)*I)/a;
m=M(i,j);
if m>0
text(V(1),V(2),num2str(m));hold on
end;
end;
end;