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

If 2 is the only base for which you want to do this, and you don't actually need the prime factorization of the remaining odd part, then there's a much faster way to do this. I mention this because factoring out the highest-possible power of 2 is a very common preliminary step in number-theory algorithms; so I'm guessing that that may be your goal. Here's how:

OddPart:= (n::posint, x::algebraic)-> 
    (k-> x^k*integerdivq2exp(n,k))(Bits:-FirstNonzeroBit(n, 'higher'))
:
#Random test case:
(p1,p2):= nextprime~(['rand(2^64..2^65)()'$2])[];
      p1, p2 := 24421176063119381659, 21314142271437361117
N:= 2^rand(9..99)()*p1*p2;
        N := 266504407575115289384393504182133531188736

CodeTools:-Usage( ifactor(N) );
memory used=14.24MiB, alloc change=0 bytes,
cpu time=172.00ms, real time=175.00ms, gc time=0ns

          9                                              
       (2)  (21314142271437361117) (24421176063119381659)

CodeTools:-Usage( OddPart(N,x) );
memory used=1.49KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

                                                    9
           520516421045147049578893562855729553103 x 

Your pde2 is not explicitly an equation because it doesn't consist of two expressions separated by =. Thus, it has neither a left side nor a right side. Thus, the command lhs (which stands for left-hand side) can't be used. This error is purely due to syntax, not math. If you intend pde2 to represent an expression implicitly equated to 0 (which is a very common and totally acceptable thing to do), then simply proceed with normal(pde2) instead of normal(lhs(pde2)), and your entire worksheet completes without error.

Regarding your tanh question: Did you intend that question to be in connection to the error? If so, the error has nothing to do with math. If you intended the tanh question to be independent of the error, I can attempt to answer it. 

The simplify with side relations (see ?simplify,siderels) command can do it:

simplify(x^2+y^2, {x-y, y+a}, [x,y]);
                                

The command to clear the counter on an automatically generated global variable such as _Z1 is `tools/genglobal`. Before each solve command, including the first, do

unprotect(_Z); `tools/genglobal`(_Z, 0, 'clear'):

However, it would be a better programming practice to "improve the program" in the manner @acer suggests. One reason that this is more robust is that if you and other people are developing the same code, they might not be aware that you've reset the counter.

`tools/genglobal` is undocumented. I figured out how to use it by reading its code via

showstat(`tools/genglobal`);

Admittedly, this requires substantially more effort than reading a help page.

The following simulation, sampling over 2 million random drawings, shows that the probability that a random drawing would have no intra-family pairings is less than 2%:

Fams:= [{$1..3}, {$4..7}, {$8..11}, {$12..15}, {16}]:
ValidDraw?:= (M::list(integer[1..16]))->
    andmap(f-> f intersect {M[[f[]]][]} = {}, Fams)
: 
V:= [
    to 2^21 do 
        M:= combinat:-randperm(16); 
        if ValidDraw?(M) then M fi 
    od
]:
nops(V);
                             37276
evalf(%/2^21);
                         0.01777458191

Using formal Statistics language, the p-value is less than 2%.

A well-known formula for combinatorial derangement says that the probability that no person is paired with themself is very close to exp(-1) ~ 36.8%. The exact value is  15549624751/42268262400, which equals exp(-1) to 14 decimal places.