The most intuitive construction for this is to have a infinitely-sized matrix A which performs
$$ V(x) \cdot A = V(a^x) \qquad \qquad \text{ where } V(x)=[1,x,x^2,x^3,...] $$
Such a matrix A is called "Carleman-matrix" (see wikipedia).
The idea is, that then
$$ V(x) \cdot A = V(a^x) = V(\exp_a(x))\\
V(a^x) \cdot A = V(a^{a^x}) = V(\exp_a^{[2]}(x)) \\
\ldots $$
or
$$ V(x) \cdot A^h = V(\exp_a^{[h]}(x)) $$
thus requiring the h'th power of the matrix $A$.
Approximations to this can be made by finite truncations of A - so let's in the following denote such finite truncations with the hat-symbol: $ \hat A$ .
It is so far not difficult to implement arbitrary precision to this if the base $a$ is in the range $\small 1 < a < \exp(\exp(-1))$ with finite matrices $\hat A$ - however, who needed such a construction for the integer height tetration?
But from here one idea is apparently straightforward: we need fractional powers of A to have fractional heights of tetration. For finite-size matrices one can define fractional powers using diagonalization and take the scalar elements of the diagonal matrix to the (fractional) h'th power. So finite truncations of A and its diagonalization gives an approximation.
Let the diagonalization of the truncated matrix $ A$ be written as
$$\hat A = \hat M \cdot \hat D \cdot \hat M \, ^{-1}$$ and for notational convenience let's write $\small \hat W$ for $ \small \hat M \, ^{-1}$
then $$ \hat A^h = \hat M \cdot \hat D \, ^h \cdot \hat W = \hat M \cdot \text{diag}( (d_{k,k})^h) \cdot \hat W $$ and -in principle-
$$ \hat V(x) \cdot \hat A \,^h \approx \hat V(\exp_a^{[h]}(x))$$
Unfortunately the truncation of $A$ spoils its Carleman-characteristic, and the best approximation is to use the second column only
$$ \hat V(x) \cdot \hat M[,1] = \hat \sigma_a(x) \\
\hat V(\hat \sigma_a(x)) \cdot \hat D \ ^h \cdot \hat W[,1] \approx \exp_a^{[h]}(x)$$
See a worked example in Pari/GP in my answer at the earlier question (your question here seems to get marked as a duplicate of that earlier one)
Remark: the last procedere mimicks in principle the method of Ernst Schröder (the $ \hat \sigma _a(x)$ represents here the "Schröder-function") and can be accurate to the -in my opinion- currently best concept for the real tetration on the real axis by H. Kneser (implemented by members of the Tetrationforum) to 10 to 12 or more digits using matrix size of 32x32 . But even size 16x16 allows reasonable approximation to say 6 or 8 digits (don't know exactly from the top of my mind)
Another related concept, immediately implementable for bases a in the mentioned interval is to use Carleman-matrices for the $\exp_a()$-function but where the power series are recentered at an attractive fixpoint. Then the Carleman-matrices become triangular and diagonalization is simpler (and the results are slightly different). Also various other concepts have been discussed but are mostly far more complicated.