Finding polynomial without constant term that commutes with $f(x)=x^3+3x$

We have $$f(x):=x^3+3x$$ and you already found $$h(x)=x^5+5x^3+5x$$ with $$h(f(x))=f(h(x))$$ but the coefficient of $x$ of $h$ is not divisible by $3$. But $$g(x)=h(f(f(x)))$$ has this property because

$$g(x)=h(f(f(x)))=f(h(f(x))=(h(f(x))^3+3(h(f(x))$$ and the lowest non vanishing power of $(h(f(x))^3$ is $x^3$ and all coefficients of $x$ in $3(h(f(x))$ are divisible by $3$.

And $$g(f(x))=h(f(f(f(x)))=f(h(f(f(x)))=f(g(x))$$

(I missed that g(x):=h(f(x)) already has the required property- It is equal to the degree 15 polynomial I calculated in the second part of this post)

The polynomial $g(x)$ is

$${{x}^{45}}+45 {{x}^{43}}+945 {{x}^{41}}+12300 {{x}^{39}} \\+111150 {{x}^{37}}+740259 {{x}^{35}}+3764565 {{x}^{33}}+14945040 {{x}^{31}}\\ +46955700 {{x}^{29}}+117679100 {{x}^{27}}+236030652 {{x}^{25}}+378658800 {{x}^{23}}\\ +483841800 {{x}^{21}}+488494125 {{x}^{19}}+384942375 {{x}^{17}}+232676280 {{x}^{15}}\\ +105306075 {{x}^{13}}+34512075 {{x}^{11}}+7811375 {{x}^{9}}+1138500 {{x}^{7}}\\ +95634 {{x}^{5}}+3795 {{x}^{3}}+45 x$$ which can be factored to

$$x\, \left( {{x}^{2}}+3\right) \, \left( {{x}^{4}}+5 {{x}^{2}}+5\right) \, \left( {{x}^{6}}+6 {{x}^{4}}+9 {{x}^{2}}+3\right) \,\\ \left( {{x}^{8}}+7 {{x}^{6}}+14 {{x}^{4}}+8 {{x}^{2}}+1\right) \, \\ \left( {{x}^{24}}+24 {{x}^{22}}+252 {{x}^{20}}+1519 {{x}^{18}}\\ +5796 {{x}^{16}}+14553 {{x}^{14}}+24206 {{x}^{12}}+26169 {{x}^{10}}\\ +17523 {{x}^{8}}+6623 {{x}^{6}}+1182 {{x}^{4}}+72 {{x}^{2}}+1\right) $$

Here is a smaller one

${{x}^{15}}+15 {{x}^{13}}+90 {{x}^{11}}+275 {{x}^{9}}+450 {{x}^{7}}+378 {{x}^{5}}+140 {{x}^{3}}+15 x$

I calculated some polynomials with Maxima. The following commute with $f$ and are the only one with nonnegative leading coefficient and a highest power less or equal $15$:

$$0$$ $$x$$ $$x^3+3x$$ $$x^5+5x^3+5x$$ $$x^7+7x^5+14x^3+7x$$ $${{x}^{9}}+9 {{x}^{7}}+27 {{x}^{5}}+30 {{x}^{3}}+9 x$$ $${{x}^{11}}+11 {{x}^{9}}+44 {{x}^{7}}+77 {{x}^{5}}+55 {{x}^{3}}+11 x$$ $${{x}^{13}}+13 {{x}^{11}}+65 {{x}^{9}}+156 {{x}^{7}}+182 {{x}^{5}}+91 {{x}^{3}}+13 x$$ $${{x}^{15}}+15 {{x}^{13}}+90 {{x}^{11}}+275 {{x}^{9}}+450 {{x}^{7}}+378 {{x}^{5}}+140 {{x}^{3}}+15 x$$

Note that

$${{x}^{9}}+9 {{x}^{7}}+27 {{x}^{5}}+30 {{x}^{3}}+9 x=f(f(x))$$ and $${{x}^{15}}+15 {{x}^{13}}+90 {{x}^{11}}+275 {{x}^{9}}+450 {{x}^{7}}+378 {{x}^{5}}+140 {{x}^{3}}+15 x = h(f(x))$$ where $$h(x)=x^5+5x^3+5$$

If we call these polynomials (excluding 0) $h_1,h_3,h_5,\ldots, h_{15}$, so the index $i$ of $h_i$ is the degree of the polynomials, the following holds:

$$h_i(x)=(x^2+2)h_{i-1}(x)-h_{i-2}(x)$$

Here is an online version of Maxima. This is how I calculated the coefficients of $h_9$ by comparing the coefficients of $(f(g(x))$ and $g(f(x))$ N:9 is the degree of the polynomial g(x). The statement starting with empty checks if all coefficients of the difference $(f(g(x))-g(f(x))$ become $0$. If the result of this statement is false, the solution is invalid. This happens if we choose $N$ and even number.

f(x):=x^3+3*x;
N:9;
g(x):=x^N+sum(a[i]*x^(i-1),i,1,N);
d:f(g(x))-g(f(x)),expand$
coeff_list:makelist(coeff(d,x,i),i,3*N,0,-1)$
eqn_list:makelist(coeff_list[i],i,2,N+1)$
var_list:makelist(a[i],i,1,N)$
solutions:solve(eqn_list,var_list);
emptyp(sublist(coeff_list,lambda([t], is(t#0)))),solutions;
expression:g(x);
expression,solutions;

Here is the same algorithm implemented as SmPy scirpt.

These polynomials are related to Chebyshev Polynomials as some commenters pointed out. In the comments alex.jordan linked to A082985 and Lubin mentions the comments of achillehui of the OPs questions, which contains references to a paper on this topic by F.J.Ritt, Permutable Rational Functions.