There is a long-standing (well over 20 years), fairly well-known, formal mechanism for doing this. This can be done for any (unevaluated) function:

`print/pochhammer`:= (a,n)-> ``(a)[n];
pochhammer(1/2, n);

Unlike the proposed methods using subs or eval, the method above only affects the prettyprinted display of the function; it doesn't change the function itself or your ability to perform further computations with it.

I think that what's happening is that due to a design oversight, the overload wrapping obscures the special status of a parameter named _self in an object method. A workaround is to use a normal procedure as a wrapper of the overload. The normal procedure explicitly passes _self to the overload, and implicitly (via _rest) passes the user's other arguments.

restart;
module person() option object;
local 
    Name::string:= "*", Age::numeric, #non-static object data

    set_info_inner::static:= overload([
        proc(self::person, name::string, $) option overload;
            self:-Name:= name;
            return
        end proc,

        proc(self::person, name::string, age::integer, $) option overload;
            self:-Name:= name;  self:-Age:= age;
            return
        end proc,

        proc() error "Wrong number or type of arguments", args[2..] end proc
    ])
;
export
    #wrapper that passes _self to overload:
    set_info::static:= _self-> set_info_inner(_self, _rest),

    Name?::static:= proc(_self, $) Name end proc,
    Age?::static:= proc(_self, $) Age end proc
;
end module
:
o:= Object(person);
o:-set_info("me"):
o:-set_info("me", 20):

o:-Name?();
o:-Age?();
Name?(o), Age?(o);
              o := Object<<person,1582091427520>>
                              "me"
                               20

                            "me", 20

o:-set_info("me", "you", 20);
Error, (in person:-set_info_inner) Wrong number or type of arguments, me, you, 20

  By the way, I haven't seen you on MaplePrimes for a long time. I'm glad that you're back.

Try adding this option to whatever plot command you used to generate that plot:

axis[1]= [tickmarks= [["0 - 2", "3 - 5", "6 - 8", "9 - 11"], rotation= Pi/2]]

If that doesn't work, please post the actual plot command you used; I'm sure that some minor adjustment will work. The font for the tick labels can also be changed.

Replace sin(x) with sin(x[i]), and replace sin*t[j] with sin(t[j]).

Here are two more ways: If you know that you want the index of the maximum element, use

max[index](R);

To get the index of any element, use

member(96, R, 'p'):  p;

where p is any variable.

Let T be any type. Then

• Type [T] matches lists with exactly 1 element, which must be type T;
• Type list(T) matches lists with any number of elements (including 0), all of which are type T.

If you want to forbid the case of empty lists, you could make the declaration

m::And({Matrix, list(Matrix)}, Not([])) 

or

m::And({thistype, list}(Matrix), Not([])).

thistype is the identity function for types (not to be confused with the identical functional type).

Assuming that the 5 terms that you call Lr represent the beginning of an infinite power series in r, then that's a simple geometric series that converges iff abs(r) < 2/Pi, and its exact value is 

2*r/(1 + (2/Pi/r)^2),

or, equivalently,

2*r*(Pi*r)^2/((Pi*r)^2 + 4).

Since this exact value is easily obtained and is a degree-[3,2] rational function, it doesn't make sense to me to use a pade approximation.

Even if you truly only want to consider those 5 terms, they can be considered as the difference of two infinite geometric series to obtain

2*r*(1 + 1/R^5)/(1 + R), where R = (2/Pi/r)^2,

which is mathematically identical to your original 5 terms, not an approximation (with removable singularities for r in {-1, 0, 1}*2*I/Pi, but no inequality restriction on abs(r)).

I often need to filter a set/list of solution sets/lists to separate those that contain trivial equations. It can be done like this:

s:= [{f0(t, r) = 0, n(t) = n(t)}, {f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}];
selectremove(t-> ormap(evalb, t), s);
     [{f0(t, r) = 0, n(t) = n(t)}],  [{f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}]

This code would be identical regardless of whether the solutions are sets or lists. 

[This Answer was based on the original Question alone, before the OP posted any code. The Replies were written after the code was posted.]

I suspect that this is the correction for your problem: Let's suppose that you have a procedure P that creates a table assigned to a variable T that is a local of P and you want that table to be the return value of P. Then due to the last name evaluation property of T, you should make the return value eval(T) rather than simply T. You should do this inside P. No second argument for the eval is necessary, but after you get this working with the simple eval(T), you might be able to make a slight efficiency improvement using an integer second argument for the eval depending on the what type of things are stored in T.
 

I suspect that this is the cause of your problem: Suppose that P returns T without using eval. You make  your call P(...), and the displayed output makes it seem as if T has been returned. But that T is still the T local to P. That's called an escaped local. (Although "escaped" might make them sound dangerous, they can be useful.) It is a distinct variable from the global T, which can use outside of P, such as in a worksheet. The two Ts have no connection whatsoever; they merely have names that are coincidentally spelled the same. I suspect that you've been trying to use something akin to T[k] from your worksheet level.
 

You can assign the returned table to any variable: If you use eval as I described in the first paragraph, and you make your call to P as T:= P(...), then you can index T at the worksheet level. But you can use any variable for this purpose; it needn't be T.
 

I am certain that these things have nothing to do with your problem: It makes absolutely no difference where the procedure P is defined (.mpl file, library archive, or directly in your worksheet). It makes no difference where you call P from (worksheet or from another procedure).

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: