2 of 2 alternate approach

PARI/GP, 4127 digits

(104127-1)/9 + 2*10515

This is a fairly straightforward search: check only prime digit lengths, then compute the possible primes to use, then iterate through all possibilities. I special-cased the common cases where there are 0 or 1 suitable prime digits to use.

supreme(lim,startAt=3)={
forprime(d=startAt,lim,
my(N=10^d\9, P=select(p->isprime(d+p),[1,2,4,6]), D, n=1);
if(#P==0, next);
if(#P==1,
for(i=0,d-1,
if (ispseudoprime(D=N+n*P[1]), print(D));
n*=10
);
next
);
D=vector(#P-1,i,P[i+1]-P[i]);
for(i=0,d-1,
forstep(k=N+n*P[1],N+n*P[#P],n*D,
if (ispseudoprime(k), print(k))
);
n*=10
)
)
};
supreme(4200, 4100)


This took 36 minutes to compute on one core of a fairly old machine. It would have no trouble finding such a prime over 5000 digits in an hour, I'm sure, but I'm also impatient.

A better solution would be to use any reasonable language to do everything but the innermost loop, then construct an abc file for primeform which is optimized for that particular sort of calculation. This should be able to push the calculation up to at least 10,000 digits.

Edit: I implemented the hybrid solution described above, but on my old machine I can't find the first term with >= 10,000 digits in less than an hour. Unless I run it on something faster I'll have to change to a less lofty target.