This is true. First, let $f(x)=x^n+5x+3$. For any polynomial $g(x)$ of degree $n$, we define its reciprocal polynomial as $\tilde{g}(x)=x^ng(\frac{1}{x})$. Thus, $\tilde{f}(x)=3x^n+5x^{n-1}+1$.
This proof follows the method described in Irreducibility of $x^n-x-1$ over $\mathbb Q$. Here we show two things:
- $f(x)$ and $\tilde{f}(x)$ do not share a root
- for all integer polynomials $k(x)$ of degree $n$ such that $f(x)\tilde{f}(x)=k(x)\tilde{k}(x)$, we have either $k(x)=\pm f(x)$ or $k(x)\pm \tilde{f}(x)$
The referenced result implies that $f(x)$ is irreducible. This method is most effective when $f(x)$ has few small coefficients (ideally $\pm 1$). Since this is not our case, proving the second point requires more extensive case analysis. It would be possible to do by hand, but it would require to check about 60 cases - so I use Python to automate this process.
$f(x)$ and $\tilde{f}(x)$ do not share a root
For the sake of contradiction, assume there exists $\alpha\in\mathbb{C}$ such that $f(\alpha)=\tilde{f}(\alpha)=0$. Then
$$
0=5f(\alpha)-\alpha(\tilde{f}(\alpha)-3f(\alpha))=15\alpha^2+33\alpha+15.
$$
Hence the irreducible polynomial $5x^2+11x+5$ divides $x^n+5x+3$ in $\mathbb{Z}[x]$, impossible since $5\nmid 1$. Hence no such root exists.
$f(x)\tilde{f}(x)=k(x)\tilde{k}(x)$ implies $k(x)=\pm f(x)$ or $k(x)=\pm \tilde{f}(x)$
Let $k(x)=b_0+b_1x+\dots+b_nx^n$ and
$$
f(x)\tilde{f}(x)=3+5x+15x^{n-1}+35x^n+15x^{n+1}+5x^{2n-1}+3x^{2n}=k(x)\tilde{k}(x).
$$
Comparing the coefficients, we get restrictions on coefficients $b_i$. Looking at $x^n$:
$$
35=b_0^2+b_1^2+\dots+b_n^2\tag{*}
$$
Without loss of generality, we can replace $k(x)$ with $\tilde{k}(x)$ and multiply both by $-1$ to achieve $b_0=3$ and $b_n=1$. So (*) gives
$$
25=b_1^2+\dots+b_{n-1}^2\tag{**}
$$
Looking at $x^1$ we get
$$
5=b_0b_{n-1}+b_1b_{n}=3b_{n-1}+b_1.
$$
Hence $b_1\equiv 2\pmod 3$ and from $(**)$ we get $|b_1|\leq 5$. Thus $b_1\in \{-4,-1,2,5\}$. To each value of $b_1$ we express $b_{n-1}$ to find
$$
(b_1,b_{n-1}) \in \{(-4,3),(-1,2),(2,1),(5,0)\}.\tag{***}
$$
Subtracting $b_1^2+b_{n-1}^2$ from $(**)$, we repeat the same process for each case and inspect the coefficients of $x^2$ to determine $(b_2,b_{n-2})$, $x^3$ to determine $(b_3,b_{n-3})$, and so on. Assuming these coefficients are zero (which they are up until the term $15x^{n-1}$), we apply the method recursively and when the condition $(**)$ reduces to $b_l^2+\dots+b_{n-l}^2=0$, we know all of the coefficients of $k(x)$.
The above process yields 58 such polynomials (or more precisely, polynomial families, since $n$ is a variable). However, when additional constraints are applied, requiring more coefficients to be zero (up to $x^{l+2}$), only a single case remains:
$$k(x)=x^n+5x+3.$$
Thus $k(x)=f(x)$, and the result follows.
Throughout this process, the maximal value of $l$ is $16$. Based on the order of exponents, we require $l+2<n-1$ and $l<n-l$. Thus our argument holds for $n>32$. The remaining polynomials can be checked individually, as you have done.
Here are more technical details of the recursive algorithm described above.
Inputs: $(b_0,b_n), (b_1,b_{n-1}),\dots, (b_l,b_{n-l})$ with $b_0=3, b_n=1$
Let $K=35-b_0^2-\dots-b_l^2-b_n^2-\dots-b_{n-l}^2$
- If $K=0$, then
$k(x)=b_0+b_1x+\dots+b_lx^l+b_{n-l}x^{n-l}+\dots+b_nx^n$ is a
candidate. For $i=1,2,3$ we check if coefficient at $x^{l+i}$ in
$k(x)\tilde{k}(x)$ is $0$, i.e. $$ b_lb_{n-i}+\dots+b_ib_{n-l}=0. $$
If this holds, we report $k(x)$.
- If $K>0$, let $B=-\sum_{i=1}^{l}b_ib_{n-l-1+i}$, and by looking at
coefficient $x^{l+1}$ we have $b_{l+1}\equiv B\pmod 3$ and also
$|b_{l+1}|\leq \lfloor\sqrt{K}\rfloor$. For each such $b_{l+1}$ we
let $b_{n-l-1}=(B-b_{l+1})/3$. We repeat the process with the
additional $(b_{l+1},b_{n-l-1})$.
Starting from the four initial cases in $(***)$, this process outputs:
[(3, 1), (5, 0)]
max l = 16
(see Python implementation)
This corresponds to $k(x)=x^n+5x+3$ and concludes the irreducibility for $n>32$.