For example, to do it for year 2024:

select(m-> Calendar:-DayOfWeek(2024, m, 13) = 6, [$1..12]);
                          [9, 12]

So the 9th and 12th months have a Friday-the-13th, i.e., September and December. The 6 represents Friday (considered the 6th day of the week).

Use 

if indets(r, suffixed('_')) <> {} then . .

You could also change '_' to _Z, etc.

There are two anomalies happening in your example:

1. Left-associativity: All infix operators beginning with a single & are left-associative (and those beginning with && are right-associative). So 2 &^ 2 &^ 2 &^ 2 is interpreted as ((2 &^ 2) &^ 2) &^ 2 rather than 2 &^ (2 &^ (2 &^ 2)). This happens in the kernel before mod sees the expression.

2.Distributivity: mod doesn't distribute over exponentiation because it's usually not mathematically valid. In other words, a &^ b mod c isn't interpreted as (a mod c) ^ (b mod c) mod c. An example for which it isn't valid is a=2, b=4, c=3. 

I wouldn't use testeq at all. It's returned FAIL every time that I've ever tried to apply it to a transcendental function. But plain simplify or is works fine here. All of these return the information that you want:

is(expr=x), evalb(simplify(expr)=x), evalb(simplify(expr-x) = 0);

Here are two more-efficient algorithms for the sum of the proper divisors. The first of these (SPD2 below) is much faster than your original; the second (SPD3 below) is much faster than that.

I think that this will all work in your Maple 13; if not, let me know.

restart
:
(* Sum of Proper Divisors (original version):
---------------------------------------------
This is your original procedure. The only changes I made were removing things
not directly related to getting the answer, such as print statements.
*)
SPD0:= proc(n)
local Sum1, b;
    Sum1:= 1;
    for b from 2 to ceil(n/2) do
        if floor(n/b) = n/b then Sum1:= Sum1 + b fi
    od;
    Sum1
end proc
:
(* Sum of Proper Divisors (version 1):
--------------------------------------
This is the same algorithm as yours, recoded. The changes that I made were 
to use operators that are more efficient for integer arithmetic: add, modp,
and iquo instead of `for`, ceil, `/`, and floor.
*)
SPD1:= proc(n::posint)
local d;
    1 + add(`if`(modp(n,d) = 0, d, 0), d= 2..iquo(n,2))
end proc
:
(* Sum of Proper Divisors (version 2):
--------------------------------------
This algorithm uses this idea: 
Let n be a positive integer. Let s = floor(sqrt(n)). For every divisor d of n, 
q = n/d is also a divisor, and d <= s iff q >= s. If s^2 = n, we need to adjust
for the double counting of s.
*)
SPD2:= proc(n::posint)
local d, q, s;
    s:= isqrt(n);  #same as floor(sqrt(n)) but faster
    #irem gives the integer remainder and quotient in one step:
    add(`if`(irem(n, d, 'q') = 0, d+q, 0), d= 2..s)  +  1 - `if`(s^2=n, s, 0)
end proc
:
(* Sum of Divisors of a Prime Power:
------------------------------------
For prime p, the divisors of p^e are clearly
    {1, p, p^2, ..., p^e}.
The sum of those is a simple geometric sum:
    sum(p^k, k= 0..e) = (p^(e+1) - 1)/(p - 1).
*) 
SDPP:= (p::prime, e::posint)-> (p^(e+1) - 1)/(p-1)
:
(* Sum of Proper Divisors (version 3):
--------------------------------------
This uses the classic formula for the sum of divisors, then subtracts
the number itself.

This uses ifactors, which returns the same information as ifactor, but in a
list of lists format:
   ifactors(n) = [-1 or +1, [[p1,e1], [p2,e2], ..., [pk,ek]]],
such that the prime factorization is 
    n = p1^e1 * p2^e2 * ... * pk^ek.
*)
SPD3:= proc(n::posint) 
local P; 
    mul(SDPP(P[]), P= ifactors(n)[2]) - n 
end proc
: 
#Test all 4 on your example:
SPD||(0..3)(945);
                       975, 975, 975, 975

#Construct a much larger random example:
seq('rand(2..99)(), nextprime(rand(10^5..10^7)())', k= 1..9);
    94, 8633311, 97, 9909901, 60, 1584343, 70, 3888067, 97, 3079871, 
    91, 5954371, 19, 8014679, 57, 682901, 84, 4397777

n:= `*`(%);
n := 7153933542604628643414930949912947237856549377795642829298491576067369227200

CodeTools:-Usage( SPD3(n) );
memory used=9.62MiB, alloc change=0 bytes, 
cpu time=94.00ms, real time=160.00ms, gc time=0ns

28515792740732679749844128611744262891259901154773757373417700467974788852800

Here's something. You'll probably want some adjustments, but it's a good start.

Pl:= plots:  PT:= plottools:
Pl:-display(
    PT:-line~(
        [[0,1,1], [0,-1,-1], [5,1,-1]], [[0,1,-1]$3], 
        linestyle= dash, color= red
    )[],
    PT:-line~(
        [
            [0,1,1], [5,-1,1], [0,-1,-1],
            [0,-1,-1], [5,-1,1], [5,1,-1],
            [0,1,1], [5,-1,1], [5,1,-1]
        ], 
        [[0,-1,1]$3, [5,-1,-1]$3, [5,1,1]$3],
        color= red
    ),
    PT:-line~(
        [[0.5,-1.3,-1.3], [3,-1.3,-1.3]], 
        [[2.5,-1.3,-1.3], [5.0,-1.3,-1.3]],
        color= blue
    ),
    Pl:-textplot3d([2.75, -1.3, -1.3, `ℓ`], color= blue),
    Pl:-textplot3d(
        [2.5, 0, 0, typeset(Q(x,z,t) = Q[0]*delta(x - v*t)*sin(p*z))],
        color= black, align= below
    ),
    axes= normal, axis[1,2,3]= [color= blue, tickmarks= 0],
    scaling= constrained, view= [0..6, (-2..2)$2],
    labels= [x,y,z], projection= 0.7, orientation= [-91,65,-15]
);

Replace your ID with thistype.

The way that you present the function f(x,y) suggests to me that you want a plot of a vector field rather than a pair of surfaces. Otherwise, my suggestion is very similar to @mmcdara :

Explore(
    plots:-fieldplot([-a*x/(1+y^2), x+b*y], (x,y)=~ -5..5),
    (a,b)=~ -2.0 .. 2.0
);

Change the argument of arctan to D(z)(x) or D(z)(0). Since x=0 at this point in your code, your derivative expression is asking for the "derivative with respect to 0", which, of course, is nonsense. Using D(z)(x) tells it to first take the symbolic derivative of z, and then evaluate that at x (which could be simply a number, like 0 in this case).

Unrelated issue: A design "wart" of Maple's units implementation is that 0 cannot "hold onto" its units, because the Maple representation is number x unit. If the number is 0, that simplifes to just 0. This may or may not cause you problems later on; it's just something to be mindful of. But using 0 with a unit as you did will not in and of itself cause you any problems, and it helps make your code more clear.

From the menus, go to Tools => Options => Interface and check the box "Remember plot attributes when re-executing", then click Apply Globally.

You can use

C||(1..numelems(V)):= seq(V);

Edit: Thanks to @nm for noticing that I'd erroneously used nops instead of numelems.

Maple's numerical ODE solver can solve BVP eigenvalue problems. It will give a solution, if it can find one. It doesn't guarantee uniqueness or minimality.

Download eigvBVP.mw

You have nested add commands of the form

add(add(add(..., r= 0..k), p= 0..r), l= 0..p)

add(add(add(add(..., r= 0..k), p= 0..r), l= 0..p), j= 0..l)

You have the nesting order completely backwards. They should be

add(add(add(...., l= 0..p), p= 0..r), r= 0..k)

add(add(add(add(..., j= 0..l), l= 0..p), p= 0..r), r= 0..k)

The regular plot command can be used for this. You just need to group the X and Y into a single argument:

plot(<X1 | Y>, style= point, symbolsize= 15, useunits= Unit~(['m', 1]));

Your only have 1 equation (Eq 21), so it can only be solved for 1 variable at a time. Your stated "solutions" have solutions for 4 variables: omega, Psi, kappa__0, and h__1. If we take 3 of these as "given", we can easily solve for the 4th. In the worksheet below, I solve for omega.

Download Xi.mw