Having a variable, a name, that is assumed to be integer is not the same as having a value that is actually an integer (i.e., has type integer). The command numtheory:-bigomega wants an actual integer.

Now let's learn about premature evaluation, which is the essential difference between

sum(numtheory[bigomega](k), k= 2..9)

and

add(numtheory[bigomega](k), k= 2..9)

The command sum is like the vast, vast majority of other Maple commands: It evaluates its arguments before they are passed. That means that the above sum command attempts to evaluate numtheory[bigomega](k) before k has been assigned any numeric value. That causes an error.

On the other hand, the commands seq, add, and mul have special evaluation rules: Their index variable is assigned values from the prescribed range before the first argument is evaluated.

To some extent, premature evaluation can be controlled with unevaluation quotes; however, those are often difficult or impossible to use. In the above case, sum('numtheory[bigomega](k)', k= 2..9) works correctly, but it's better to just use add.

When you get an error message that says that a certain argument is expected to be a certain type, that always means that that argument must actually be that type; that error can never be corrected by making assumptions.

Your definition of N is correct. It's obvious that the midrange is approximately 10 and that the supposed answer is way off.

To get the min and max of a 2D plot, first assign the plot to a variable:

P:= plot(N, 0..20);

then

(min,max)(op([1,1], P)[.., 2]);

      9.49146814577784781, 11.3050481854892997

You can simply do

1 +~ ListTools:-PartialSums([1, 3, 4, 50, 10]);

You can't use the same variable, i, as both the index for the for loop and the index for the sum command. Simply use a different variable.

There are several other improvements that you could make, which I'll detail in a Reply.

Do

Matrix((2,2), (i,j)-> eval(p[j], solution[i]));

a) If you have actual data points, you should use Statistics:-PowerFit. It'll do the log transformation and find the closest exponent for you.

b) If R = a*P^(2/3), then R/P = a*P^(-1/3). From there, the two limits are trivial.

c) You have a lot of questions. Are you trying to do your homework for the entire semester in one night? You should indicate/show that you understand the answers that I already gave before proceeding.

Like this

restart:
K:= 100:  N__0:= 20:
NP:= R-> rsolve({N(t+1) = R*N(t)/(1+(R-1/K)*N(t)), N(0)= N__0}, N(t), makeproc):
plot(
   [seq(['[k, NP(R)(k)]' $ k= 1..10], R= [2, 5, 10])],
   legend= (R=~ [2, 5, 10]),
   labels= [t,N]
);

If it appears as a straight line on a log-log plot, then the equation must be as I've given as Meq, for some yet-to-be-determined values of k and C. The two given data points are enough to solve for these.

restart:
Meq:= t = exp(k*ln(T)+C):
S:= solve({eval(Meq, [t=3, T=20]), eval(Meq, [t=20, T=5])});  
plots:-loglogplot(eval(t, eval(Meq, S)), T= 1..30);

Do

restart:
R:= k*(a-x)*(b-2*x)^2:
a:= 5:  b:= 6:  k:= 0.3:
plot(R, x= 0..min(a,b/2), labels= [x,'R']);

You can also plot x as a function of time like this:

Sol:= dsolve({diff(x(t),t) = subs(x= x(t), R), x(0)=0}, numeric):
plots:-odeplot(Sol, [t,x(t)], t= 0..min(a,b/2));

In Maple's plaintext input (aka 1D input or Maple Input), the juxtaposition of parenthesized expressions does not imply their multiplication. In other words (a-1)(a-2) doesn't mean (a-1)*(a-2). Unfortunately, it does mean something, so it's not a syntax error. What it means is rather obscure and not worth discussing here.

Do you really need a Groebner basis? Or do you just need a solution? The equations can be easily solved by solve. Here, I just use equations 17-44, because I didn't feel like correcting the first 16, and it's trivial to see that the given solution satisfies the first 16.

solve({e||(17..44)});
      {a = 2, b = 3, c = 4, d = 1, e = 1, f = 4, g = 3, h = 2, i = 3, j = 2, k = 1, l = 4, m = 4, n = 1, o = 2, p = 3}

Here's a pair of procedures for it. They'll work in any base, although the results will always be reported in base 10. The second parameter, n, of KaprekarSeq is the maximum number of terms of the sequence to report. This is needed because it's not guaranteed that the sequence will converge to a fixed point. It's likely that it will, but not guaranteed.

KaprekarSeq:= proc(x::posint, n::posint, b::posint:= 10)
local j, r, kx;
   r[0]:= x;
   for j to n-1 do 
      kx:= Kaprekar(r[j-1], b);
      if kx = r[j-1] then break end if;
      r[j]:= kx
   end do;
   convert(r, list)
end proc:

Kaprekar:= proc(x::posint, b::posint:= 10)
local X:= convert(x, base, b), N:= b^nops(X), S:= sort(X);
   convert(S, base, b, N)[] - convert(ListTools:-Reverse(S), base, b, N)[]
end proc:

Example usage:

KaprekarSeq(1234, 9);

       [1234, 3087, 8352, 6174]

I don't know how to solve them all together, or if it's possible at all, but we can easily use dsolve to solve (1) & (2) to get psi__1 and psi__2 as functions of x and likewise with (3) & (4) to get functions of t. Like this:

eq1:= diff(psi__1(x), x) = -I/2/lambda*(rho*psi__1(x) + alpha*exp(I*(x+rho/2*t))*psi__2(x)):
eq2:= diff(psi__2(x), x) = -I/2/lambda*(-rho*psi__2(x) + alpha*exp(-I*(x+rho/2*t))*psi__1(x)):
Sol__x:= dsolve({eq1, eq2});

eq3:= diff(psi__1(t), t) = I*lambda/2*psi__1(t) - alpha/2*exp(I*(x+rho/2*t))*psi__2(t):
eq4:= diff(psi__2(t), t) = -I*lambda/2*psi__2(t) + alpha/2*exp(-I*(x+rho/2*t))*psi__1(t):
Sol__t:= dsolve({eq3, eq4});

Do

P:= sort(R, output= permutation);
L[P]; M[P];

One option is to do this:

plots:-contourplot3d(
   f(x,y), x= 0..1, y= 0..1, contours= [0.1],
   thickness= 5, color= red
);

Why do you treat u and all of its derivatives as independent variables? The only independent variables are x, y, t. You can just do this:

a:= -4*D[1,1](u) + _F1:
convert(D[1](a)(x,y,t), diff);

This produces your expected output.