When I read @dharr 's proposed GAMMA workaround this morning, before @acer looked into it, I knew that it couldn't possibly work in the desired way (i.e., under evalhf) because GAMMA(401.) is still well above the maximum value of a hardware float. But replacing the three factorials with a single call to binomial does work under evalhf:

F:= (n,k,p)-> 100*binomial(n,k)*p^k*(1-p)^(n-k):
evalhf(F(400, 5.1, 1/100));
                      15.1978021802672245

CodeTools:-Usage(plot(F(400, q, 1/100), q= 0..20, numpoints= 1000));
memory used=2.88MiB, alloc change=0 bytes, cpu time=31.00ms, real time=40.00ms, gc time=0ns

int(diff(y(x), x), x= 1..2, continuous)

If you use option tryhard, then it uses a completely different algorithm and that limitation doesn't apply.

I have a completely different approach that doesn't use any strings or character-based representations, works in any base, works for any desired "digit" count, and is faster than the fastest procedure from the original Question by a factor of 2 (although it doesn't quite make it to the desired "quarter of a minute").

Here's the procedure:

(* 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 sequences of
	base-b "digits" with the desired "digit" 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 cther charcater-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::nonnegint:= 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), repunit:= 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:=
		if D1 = 1 then ()-> bs minus {d}
		else ()-> map(`minus`, [`if`(C[1]=0, d0, bs), mids, `if`(C[D]=ds, d1, bs)], {d})[]
		fi,
	Isprime:= eval(`if`(p1 < 5*10^9, gmp_isprime, isprime))  #*1 
;
	# Short-circuit returns for corner cases:
	if D1 < 0 then return {} fi;
	if D1 = 0 then
		return 
			if D>1 then `if`(p0 <= repunit and repunit <= p1 and Isprime(repunit), {repunit}, {}) 
			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:= repunit - add((BC:= [seq](B[[seq](C)])));
			(
				for d in bs minus `if`(C[D] = ds, {}, {0}) do
					Rep:= d*rep;				
					(
						for k in It:-CartesianProduct(FactorSelector()) do               
							p:= Rep + inner([seq](k), BC);  #*2
							if p0 <= p and p <= p1 and Isprime(p) then p else fi
						od
					)
				od
			)
		od
	}
end proc

(* Footnotes:
	*1: gmp_isprime is an undocumented builtin procedure that is faster 
	than isprime for most small cases.

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

And the worksheet:

Download DigitCountPrimes.mw

I am vehemently opposed to any so-called "mathematics" based on the base-10 digits of numbers without generalizing it to all bases. To me, it's "junk math". It looks like OEIS will admit many such worthless sequences.

Your syntax is already perfectly correct, including your placement of derivatives in boundary conditions. To finish solving, replace degree with Pi/180, and set a specific numeric approximation for infty in the Parameters (I used 5). Then do

dsolve(eval({eq1, eq2, eq3, ics, bcs}, {Parameters}), numeric)

Your system of ODEs is very simple: four linear equations with constant coefficients for the homogenous parts and with the forcing functions all of the form C[k]*exp(-t) for 4 constants C[k]. The solutions are all very smooth, monotonic functions. Thus, there's very little error even for rk1 (a.k.a. Euler's method).

These are your equations:

Your can simplify your procedure substantially while correcting the error that you asked about. In the following, I use the optional argument J to list the indices that you want to perform this operation on (because it's not clear to me where your index-pattern 1, 2, 4, 5, 7 comes from). The index -1 refers to the last vector position regardless of the actual number of entries.

oneStep_egcd2:= proc(prev::Vector, curr::Vector, J::list(posint):= [1,2,4,5,7])
    if curr[-1] <> 0 then
        (prev[], curr[J]):= (curr[], prev[J] - curr[J]*iquo(prev[-1], curr[-1]))       
    end if
end proc
:
oneStep_egcd2(prev, curr);

Maple's multi-assignments happen essentially simultaneously without needing a temporary intermediary variable. For example, compare these:

(a,b):= (1,2):  (a,b):= (b,a):  a, b;  #multi-assignment
                              2, 1

(a,b):= (1,2):  a:= b: b:= a:  a, b;   #sequential assignments
                              2, 2

I think that I have a better algorithm for this. I think it has a much better average-case time. It guarantees return of the minimal possible y (if there's any y larger than its corresponding x). It does this by stepping sequentially through y, stopping at the first that's greater than its corresponding x. And, it's a much simpler algorithm theoretically.

MinY:= proc(Coeffs::list(rational), m::list(posint))
local c, M:= ilcm(m), g, X:= 0, Y:= 0;
    `mod`:= modp; #Use nonnegative remainders only
    M/= (g:= igcd((c:= chrem(Coeffs, m)), M)); #chrem is Chinese Remainder Algorithm
    c:= g/c mod M; 
    to M do until (X:= X+c mod M) < (Y+= g);
    `if`(X<Y, [X,Y], FAIL)
end proc
:    
# -x_coeff/y_coeff for each equation:
Coef:= [-154/69, -13/716, -23/3059, -2295/4522, -6479/5396]:

# Moduli: The 3rd and 4th have their power reduced because the coeffs are reducible.
M:= [7^3, 13^3, 23^2, 17^2, 29^3]
:
seq(CodeTools:-Usage(MinY(Coef[..k], M[..k])), k= 1..nops(M)); 
memory used=2.80KiB, alloc change=0 bytes, 
cpu time=0ns, real time=2.00ms, gc time=0ns

memory used=3.09KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=12.67KiB, alloc change=0 bytes, 
cpu time=0ns, real time=0ns, gc time=0ns

memory used=16.86KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=95.96MiB, alloc change=0 bytes, 
+cpu time=2.55s, real time=2.03s, gc time=796.88ms

      [7, 49], [691, 819], [15893, 18837], [1103, 26390], [186957261, 190687133]

Array indices must be integers. You're trying to use Pi/3 as an index.

A fairly simple power series for the incomplete upper Gamma function is

GAMMA(a,z) = GAMMA(a) - z^a*Sum((-1)^n*z^n/n!/(n+a), n= 0..infinity)

provided that a is not a nonpositive integer. Although this series involves a possibly noninteger power of z, that can be factored out of the power series proper, as shown.

Here's my alternative version of MemoryInUse that runs 68 times faster on my computer:

MemUse:= proc() local r; add(r[3], r= kernelopts('memusage')) end proc:

This measures exactly the same thing as MmaTranslator:-Mma:-MemoryInUse. There is something akin to a "Heisenberg Uncertainity Principle" involved here: It's impossible to measure memory usage without having some effect on it. And of course that's compounded by the fact that the device doing the measuring is the same as the device being measured. So don't expect identical numbers.

It's not incorrect for the memory usage number to sometimes decrease. This is due to "garbage collection"---the recycling of memory addresses whose contents the kernel determines will no longer be needed. Every time that you see a change in the "Memory: " or "Time: " number on the "status bar" (bottom border of your Maple window), a garbage collection has occurred.

Although there is some possibility that "remember tables" (which are the things that the forget command forgets) are causing your memory trouble, I think that that chance is only slight. So you should focus on other causes for now.

I think that you've gotten the mistaken idea that remember tables are primarily used with trig functions. Actually, any piece of Maple code could involve hundreds of remember tables, most of which will be associated with procedures that the user didn't directly call. But, most such remember tables are cleared at every garbage collection (this is called option system).

Please attach your code. If you're having trouble doing that, please describe exactly what happens.

The Maple language doesn't have a builtin operator for logical equivalence. I guess that the 2D Input parser translates a variety arrow-type logical symbols into the language's implies operator (which prettyprints as a character like =>).

But, any character that has an HTML code (and thus begins with `&...` in ASCII) can be used as an inert infix operator, so what you want prettyprinted can be achieved by

`⇔`(not (P and Q), ``(not P) or ``(not Q));

The quotes on ``(not P) and ``(not Q) are needed to prevent automatic simplification.

You may want to get rid of some parentheses from the above (at which point the expression will only be good for display, not computation). That can be done like this

use T=Typesetting in
    T:-mrow(
        T:-Typeset(not (P and Q)), 
        T:-mo(`⇔`), 
        T:-Typeset(``(not P) or ``(not Q))
    )
end use;

The problem initially seems nonlinear, but it can be easily converted (completely) to a linear-regression problem. By "completely", I mean that the following is the complete original problem, not just some linear approximation of it. The linear form gives you much more control.

Download SolveBlockQuestion.mw

Look closely at your solve command:

slove([a11*x+a12*y+a13*z=b1,a21*x+a22*y+a23*z=b2,a31*x+a32*y+a33*z=b3[,[x,y,z]);

Surely you can figure out what's wrong with it.

In addition to an `evalf/constant/...` procedure, there are a few more details that you should add to a formal symbolic constant's definition to help it work with the type system and the assume/property/is system.

restart:
`evalf/constant/K`:= ()-> evalf(sqrt(2)*sin(sqrt(17))):
addproperty(K, {RealRange(Open(-2), Open(-1))}, {BottomProp});
constants,= K;  #Or constants:= constants, K
protect(K);