Although this doesn't answer your question, it might still be helpful:
Using Mathematica (and an exhaustive amount of computing power), I have checked every number $(9n)!+n!+1$ for $n\le 2000$ with no prime found.
Finding a prime now seems very hard. For example, if we were to model a number $n$ „being prime“ as a Bernoulli random variable with parameter $\frac{1}{\ln(n)}$ [albeit motivated by the prime number theorem, this is a very rough model, for instance one could do much better already by distinguishing even and odd numbers], and if we assume the Bernoullis to be independent for $2001\le n\le 3000$ [once again very rough], then the probability of at least one success in our model is $$1-\prod_{n=2001}^{3000}\left(1-\frac{1}{\ln\left(\left(9n\right)!+n!+1\right)}\right)\approx0.00499232.$$
If you nonetheless want to continue the search for primes, here is my Mathematica code (just replace STARTHERE and STOPHERE by the lower and upper bounds of $n$ to check):
SetSharedVariable[primes, checked]; primes = {}; checked = {};
Monitor[
ParallelDo[
If[! PrimeQ[n + 1],
If[PrimeQ[(9 n)! + n! + 1], AppendTo[primes, n]]
];
AppendTo[checked, n],
{n, STARTHERE, STOPHERE}, Method -> "FinestGrained"
],
{Sort[checked], primes}
]
EDIT: I have updated the source code because we can skip numbers $n$ for which $n+1$ is a prime number as pointed out in the comments by Sil.
EDIT 2: Here is Python source code (sadly, the prime checking function of SymPy seems to be about ten times slower than that of Mathematica)
from multiprocessing import Pool
from os import cpu_count
from sympy.ntheory import primetest
from math import factorial
import time
START = 200
STOP = 200
def check(n):
num = factorial(9*n)+factorial(n)+1
if primetest.isprime(num):
print("Found prime for", n)
return n
if __name__ == '__main__':
start_time = time.time()
with Pool(cpu_count()) as p:
primes = p.map(check, list(range(START, STOP+1)))
primes = [prime for prime in primes if prime]
print("--- {} seconds ---".format(time.time() - start_time))
Best Answer
This is not a proof, but does not conveniently fit into a comment.
I'll take into account that $n=4k$ is required, otherwise $38^n+31$ will be divisible by $3$ or $5$ as pointed in the comments.
Now, if we treat the primes as "pseudorandom" in the sense that any large number $n$ has a likelihood $1/\ln(n)$ of being prime (which is the prime number density for large $n$), the expected number of primes for $n=4,8,\ldots,4N$ will increase with $N$ as $$ \sum_{k=1}^N\frac{1}{\ln(38^{4k}+31)} \approx\frac{\ln N+\gamma}{4\ln 38} \text{ where }\gamma=0.57721566\ldots $$ and for the expected number of primes to exceed 1, you'll need $N$ in the order of 1,200,000.
Of course, you could get lucky and find it at much lower $n$, but a priori I don't see any particular reason why it should be...or shouldn't.
Basically, in general for numbers $a^n+b$, the first prime will usually come fairly early, otherwise often very late (or not at all if $a$ and $b$ have a common factor).
Of course, this argument depends on assuming "pseudorandom" behaviour of the primes, and so cannot be turned into a formal proof. However, it might perhaps be possible to say something about the distribution of $n$ values giving the first prime over different pairs $(a,b)$.