#PARI/GP, 4127 digits

(10<sup>4127</sup>-1)/9 + 2*10<sup>515</sup>

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 <a href="http://openpfgw.sourceforge.net/">primeform</a> which is optimized for that particular sort of calculation. This should be able to push the calculation up to at least 10,000 digits.