@Carl Love I streamlined the above code down to 26 lines and added over 50 lines of comments. Let me know if there's any question about how this code works or about why I chose any particular technique.

(* A161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that hav exactly D
copies of any "digit" in their base-b representation. Henceforth "digit" always means base-b digit.

Defaults based on  _OEIS_ "A161786":
	The base b defaults to 10, but can be any nonunit positive integer.
	The critical digit count D defaults to 4, but can be any positive integer.

My method is
	1. Use simple combinatorial tools from the Iterator package to construct a sequence of
	digits with the desired count.
	2. Use simple arithmetic to compose the sequence into an integer.
	3. If the integer is in the desired range and prime, save it for return.

This code doesn't use strings or any other characater-based representations in any way.
It doesn't break down integers into digits; rather, it builds large integers from digits.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::posint:= 4, $)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local  #*1 
	(p0,p1):= op(ithprime~([m, m+n-1])), ds:= ilog[b](p1), D1:= ds+1-D, bs:= {$0..b-1},
	B:= Array(0..ds, k-> b^k), rep1:= add(B), C, BC, rep, d, Rep, k, n1,
	d0:= select(x-> igcd(x,b)=1, bs), mids:= bs$D1-2, d1:= {$map(iquo, p0..p1, B[ds])},
	Prime?:= eval(`if`(p1 < 5*10^9, gmp_isprime, prime)),
	Accept:= p-> if p0 <= p and p <= p1 and Prime?(p) then p fi,
	DigitSets:= rcurry(op@eval, d= ()) @ 
		`if`(D1=1, ()-> [`if`(C[1]=0, d0, bs)], ()-> [`if`(C[1]=0, d0, bs), mids, `if`(C[D1]=ds, d1, bs)])
;
	if D1 <= 0 then `if`(D1 < 0, {}, Accept~(`if`(D>1, {rep1}, bs)))
	elif 0 in d1 then n1:= NT:-pi(B[ds]); thisproc~([m, n1+1], [n1-m, m+n-n1], b, D)
	else  # main loop:
		{for C in It:-Combination(ds+1, D1) do
			rep:= rep1 - add((BC:= [seq](B[[seq](C)])));
			(for d in `if`(C[1]=0, bs, d0) minus `if`(C[D1] = ds, {}, {0}) do
				Rep:= d*rep;
				(for k in It:-CartesianProduct(DigitSets()) do
					Accept(Rep + inner([seq](k), BC))
				od)
			od)
		od}
	fi[]
end proc

(* Code comments:

Locals:
	Static:
		p0:	smallest acceptable prime;
		p1:	largest acceptable prime;
		ds:	digit count of p1, counting the units digits as "zeroeth";
		D1:	number of digit positions NOT being specially counted;
		bs:	set of possible single digits;
		B:	array of powers of b such that B[k] = b^k;
		rep1: integer with desired number of digits, all of which are 1;
		d0:	set of allowed units digits for a multidigits prime;
		d1:	set of allowed leading digits (possibly including 0) in p1..p2;
		mids: middle digits: sequence of max(D1-2, 0) copies of bs;
		
	Variable:
		C:	iterates over all D1-combinations of digit positions (0-relative)
			(C is returned by Iterator:-Combination as a sorted vector);
		BC:	sorted list of powers of b corresponding to the postions in C;
		rep: integer with D 1-digits and D1 0-digits with the positions of the 0s given by C;
		d:	iterates over the allowed values for the specially counted digits, taking into
			account whether their positions are first and/or last;
		Rep:	d*rep, thus the digits become all D d's and D1 0s in some order;
		k:	iterates over all digit sequences that can replace the 0s in Rep
			(k is returned by Iterator:-CartesianProduct as a vector in digit-position order);
		n1:	new value of n potentially used for recursive call;

	Procedures:
		Prime?:
			Replaces isprime with the more-efficient undocumented builtin gmp_isprime if all
			its inputs will be < 5 billion;
		Accept:
			If its input is in the desired range and prime, return it;
		DigitSets:
			Constructs the sequence of allowed-digit sets to pass to Iterator:-CartesianProduct.
			C and d are unassigned at the time DigitSets is assigned, but they'll be assigned
			when it's called.
			
Body:
	The 1st branch of the if-statement is to handle cases where there are no digits other than those
	being specially counted. It's easier to handle these as separate cases.

	The 2nd branch handles the case where there are fewer digits in p0 than in p1. Since the main loop
	requires that their lengths be equal, this branch uses recursion to split the cases by length.

	Main loop:
		Consider all possible D1-combinations C of digit positions:
			BC is the powers of b corresponding to those positions;
			rep has all digits 1 or 0.
			Consider all possible digit values d for the specially counted value:
				Rep is rep with 1s replaced by d;
				Consider all possible digit combinations k to replace the 0s in Rep:
					`inner` is an undocumented builtin command to compute the inner- (aka dot-)
					product of lists as if they were vectors;
					so Rep + inner(..) is Rep with its 0s replaced by the digits in k. 
*)
:		

Download DigitCountPrimes.mw

@Jamie128 Yes, it can be animated by just making a few small changes to the code above:

R:= 2:  #max r of featured surface 
spokes:= 24:  #number of radial gridlines in base

plots:-animate(
    plot3d, #the feature
    [
        [r, theta, (1+k)*(csc(theta)/r)^2], r= 0..R, theta= -Pi..Pi,
        coords= cylindrical, style= surface, shading= zhue, transparency= 0.15
    ],
    k= -1..1,
    background= [
        plot3d(  #the base
            [r, theta, 0], r= 0..3/2*R, theta= -Pi..Pi, coords= cylindrical,
            grid= [7, 2*spokes+1], color= COLOR(RGB, 0.97$3), thickness= 2.5, glossiness= 0
        ),
        plots:-textplot3d(  #the spoke labels
            [seq](
                [7/4*R, i, 0, i],
                i= (2*Pi/spokes)*~([$0..spokes-1] -~ iquo(spokes,2))
            ), 
            coords= cylindrical, font= ["times", 10]
        )
    ],
    view= [map(`*`, -1..1, 2*R)$2, 0..20], axes= framed, 
    orientation= [40,40,0], labels= [r$2, z], labelfont= ["times", 16], 
    axis[1,2]= [tickmarks= [-R, R]]
);

We need to see the code that you entered that produced that error. You can use the green uparrow on the toolbar of the MaplePrimes editor to upload a worksheet.

@sursumCorda Okay, that's a reasonable reason.

@sursumCorda I'm curious why you usually do module references as A[':-B'] rather than A:-B?

@sursumCorda 

[Edit: This comment was directed at your file-deletion version of editing one's libraries.]

You're right, and there's an easier way to confirm it: Just edit the global sequence libname:

restart:
sav_libname:= libname;
  sav_libname := 
    "C:\Users\carlj\maple\toolbox\2023\Physics Updates\lib", 
    "C:\Program Files\Maple 2023\lib", 
    "C:\Users\carlj\maple\toolbox\DirectSearch\lib", 
    "C:\Users\carlj\maple\toolbox\Glyph2\lib", 
    "C:\Users\carlj\maple\toolbox\OrthogonalExpansions\lib"

libname:= select(L-> SearchText("Physics", L) = 0, [libname])[]:
inttrans:-invlaplace(1/(s+a), s, t);
                           exp(-a t)

libname:= sav_libname:
forget(inttrans:-invlaplace);
inttrans:-invlaplace(1/(s+a), s, t);
                            exp(a t)

@acer Based on my seat-of-the-pants attempt to intepret that Matlab code, the OP might have this in mind, a slight variation of your densityplot:

plots:-densityplot(
    argument(-exp(1/(x+I*y))), (x,y)=~ -0.5..0.5, grid= [401$2],
    axes= none, colorbar= false, scaling= constrained,               
    colorscheme= ["zgradient", ["Red","Violet"], colorspace= "HSV"]
);

@Preben Alsholm I have the bug, and I confirmed that I have the latest Physics Update:

Physics:-Version();
The "Physics Updates" version in the MapleCloud is 1561 and is 
the same as the version installed in this computer, created
2023, October 20, 2:37 hours Pacific Time.

inttrans:-invlaplace(1/(s+a), s, t);
                            exp(a t)

This is probably unrelated, but I'll note that the given time and date is nearly 3 hours in the future.

@sursumCorda I decided that the D=0 case (in other words, Is there any digit that doesn't appear?) was too distinct from the the D > 0 case to handle with this same code, so I simply forbade it.

By making some minor improvements, I got the time to under a quarter minute, indeed, under 12 seconds.

(* A161786__2
This procedure returns the primes from ithprime(m) to ithprime(m+n-1), inclusive, that have
exactly D copies of any "digit" in their base-b representation.
The base b defaults to 10 but can be any nonunit positive integer.
The critical digit count defaults to 4.

My method is
	1. Use simple combinatorial tools from the Iterator package to construct a sequence of
	base-b digits with the desired count.
	2. Use simple arithmetic to compose the sequence into an integer.
	3. If the integer is in the desired range and prime, save it.

This code doesn't use strings or any other characater-based representations in any way.
It doesn't break down integers into digits; rather, it builds large integers from digits.

Comments of the form #*n refer to footnotes ar the end of the procedure.
*)
A161786__2:= proc(m::posint, n::posint, b::And(posint, Not(1)):= 10, D::posint:= 4, $)
option `Author: Carl Love <carl.j.love@gmail.com> 2023-Oct-13`;
uses It= Iterator, NT= NumberTheory;
local 
	(p0,p1):= op(ithprime~([m, m+n-1])), ds:= ilog[b](p1), D1:= ds+1-D, bs:= {$0..b-1},
	B:= Array(0..ds, k-> b^k), rep1:= add(B), rep, Rep, BC, k, d, C, p, mids:= bs$D1-2,
	d0:= select(x-> igcd(x,b)=1, bs), d1:= {$map(iquo, p0..p1, B[ds])},
	FactorSelector:=  #*1
		if D1 = 1 then ()-> `if`(C[1]=0, d0, bs) minus {d}
		else ()-> map(`minus`, [`if`(C[1]=0, d0, bs), mids, `if`(C[D1]=ds, d1, bs)], {d})[]
		fi,
	Isprime:= eval(`if`(p1 < 5*10^9, gmp_isprime, isprime))  #*2 
;
	# Short-circuit returns for corner cases:
	if D1 <= 0 then
		return
		     if D1 < 0 then {}
	    	 	elif D>1 then `if`(p0 <= rep1 and rep1 <= p1 and Isprime(rep1), {rep1}, {})
			else select(p-> p0 <= p and p <= p1 and Isprime(p), bs)
			fi
	fi; 
	if p0 < B[ds] then 
		return thisproc(m, (p:= NT:-pi(B[ds]))-m, b, D) union thisproc(p+1, n+m-p, b, D)
	fi;

	# Main code:
	{
		for C in It:-Combination(ds+1, D1) do
			rep:= rep1 - add((BC:= [seq](B[[seq](C)])));
			(
				for d in `if`(C[1]<>0, d0, bs) minus `if`(C[D1] = ds, {}, {0}) do
					Rep:= d*rep;				
					(
						for k in It:-CartesianProduct(FactorSelector()) do   
							p:= Rep + inner([seq](k), BC);  #*3
							if p0 <= p and p <= p1 and Isprime(p) then p fi
						od
					)
				od
			)
		od
	}
end proc

(* Footnotes:
	*1: FactorSelector builds the factor sets of digits for the CartesianProduct Iterators.
	Variables C and d are not defined at the time this procedure is decalred, but they'll be defined
	before it's called.
	
	*2: gmp_isprime is an undocumented builtin procedure that is faster 
	than isprime for most small cases.

	*3: 'inner' is an undocumented builtin procedure that computes the
	inner or "dot" product of 2 lists as if they were vectors.
*)
:

CodeTools:-Usage(A161786__2(10,10,10,1));
memory used=462.70KiB, alloc change=0 bytes, cpu time=16.00ms, real time=6.00ms, gc time=0ns

            {29, 31, 37, 41, 43, 47, 53, 59, 61, 67}

P2:= CodeTools:-Usage(A161786__2(10^6, 10^6)): 
memory used=1.31GiB, alloc change=0 bytes, cpu time=12.64s, real time=11.44s, gc time=1.88s

nops(P2);
                             42720

Download worksheet: DigitCountPrimes.mw

@mmcdara For the 6-digit numbers, you've selected those that have 5 occurences of a digit rather than the requested "exactly four" occurences.

@Scot Gould The command prev[]:= curr[] copies the elements of curr to the already-existing rtable represented by prev. If prev is not already an rtable, then the command does something else (very likely undesirable). The meaning of the [] when appended to an rtable is "all indices in all dimensions". The command could be replaced by prev[]:= rtable(curr).  (rtable could be replaced by Vector or copy, but these are relatively inefficient.) On the left side of the assignment, the [] is essential, just like [1] would be essential to assign to only the first element of prev.

@dharr Using null indices also works for copying:

prev[]:= curr[];

This works for any rtables, regardless of the number of dimensions (as long as it's the same on both sides) and regardless of the starting indices.

@SitecoreQedge This AI-generated Answer is much too superficial (although mostly correct), but it does have an outright error towards the end.

@Traruh Synred For me to do anything to help you with this, you need to post a complete and executed worksheet showing the occurrence of the error, and containing all the code needed to execute the worksheet.

Was it truly your goal to combine your original three 3-vectors into a single 9-vector, as shown? Or did you intend to combine them into a 3x3 matrix? Although there's nothing incorrect about the single-vector goal, the matrix goal seems a more-normal option, and it would also simplify the syntax needed to retrieve the original vectors.