When  x = 0 or 1 mod n, or x = -1 mod n and k is odd, you'll always have x^k = x mod n.  In fact x^k = x might have as many as k^2 solutions mod n.  So there are a few messages that might be easy to decrypt.  This is not usually a problem, because you're unlikely to want to send those messages. 

Yes, you can solve systems (linear or not) of several equations.
For example:

> solve({x + y + z = 3, x + 2*y - z = 2, 2*x + y + 2*z = 1});

Here's an approximation.

with(plots): with(plottools):
  p1:= rectangle([0,0],[3/2,1/3],colour=green),
          rectangle([0,1/3],[3/2,2/3],colour=red),
          rectangle([0,2/3],[3/2,1],colour=blue):
  p2:= polygon([seq([0.85 + (1/24+1/72*(-1)^j)*sin(2*Pi*j/16),
      0.5+(1/24+1/72*(-1)^j)*cos(2*Pi*j/16)],j=0..15)],colour=white):
  alpha:= arccos(32*(1/4^2 + 1/16^2 - 1/5^2)): dh1:= (Pi-alpha)/20:
  beta:= Pi - arccos(40*(1/16^2 + 1/5^2 - 1/4^2)): dh2:= (Pi-beta)/20:
  p3:= polygon([seq([3/4-1/32 + 1/8*cos(Pi+j*dh1),1/2+1/8*sin(Pi+j*dh1)],j=-20 .. 20),
              seq([3/4 + 1/10*cos(Pi-j*dh2),1/2+1/10*sin(Pi-j*dh2)],
              j=-20 .. 20)],colour=white):
  display(p3,p2,p1,scaling=constrained,axes=none);

Please don't ask me to draw the flag of Sri Lanka next.

As the error message says, pdsolve(..., numeric) only works with two independent variables.  You have three: x, y and t.

select(type, A, integer);

Note that because of the division by r, your equation is singular at r=0.  Therefore a boundary condition at r=0 may not work like boundary conditions at non-singular points.  For some equations, all solutions will satisfy the condition; for others, no solutions will satisfy it.  In this case, as Doug notes, BesselY(0,r) blows up at 0.   On the other hand, BesselJ(0,r) is finite at 0 and has first derivative 0 there, so your boundary condition at r=0 is satisfied if _C1 = 0.   Now IC will force a value (nonzero, if alpha is) for _C2.  But then the derivative at R[2] will not  be 0 unless sqrt(-rho/Dif)*R[2] is one of the zeros of BesselJ(1,z).  And there won't be any real zeros except 0 if rho and Dif are positive.

All this is a bit too subtle for dsolve.

It's unlikely that there could be a closed-form answer with n as a symbolic variable.   On the other hand, for particular integer values of n, you might be able to get an answer as a function of a (probably involving RootOf some polynomials).  You may need to make assumptions on a to ensure the antiderivative is continuous on the interval of integration.
Or if you are willing to take your chances, use the continuous option for int.

What's the question?

I think the problem is that you are assigning a value to t[q] and then using t as the independent variable in your systems of differential equations.  You should try replacing t[q] by something else that does not interfere with the variable t, e.g. you might call it tq.

What do you mean by "compute the Big oh"?   is not a particular function, it can stand for anything whose absolute value is bounded by a constant times for x near 0 (or near infinity, depending on the context). 

Your equation involves .  You could have either or , and Maple has no way of knowing which of these to use unless you tell it by writing the differential equation in one of those two forms.  

Another way:

> with(plots): with(plottools):
   a := plot([[-2, -2], [2, -2], [2, 2], [-2, 2], [-2, -2]]):
   b := plot([-1+.2*sin(t), .2*cos(t), t = 0 .. 2*Pi]):
  animate(rotate, [b, t], t = 0 .. 2*Pi, background=a);

select won't work well with a Vector, i.e. it won't remove the entries that the function returns false for, it just replaces them with NULL, so you end up with a Vector of the same size as the original, with some entries replaced by NULL.  I'd recommend converting the Vector to a list before using select.

You need to specify a value for N.  And then you have too many boundary conditions for K=0.
For K <> 0, you have a problem with the boundary condition f(0)=0 because this makes the problem singular there.  There may be solutions with such boundary conditions, but it's rather delicate - possibly too delicate for numerical methods to work well.

Where did the problem come from?  Do you have reason to believe that a solution exists?

> plots[pointplot3d]([seq(seq(seq(
    [r*sin(i*Pi)*cos(j*Pi), r*sin(i*Pi)*sin(j*Pi), r*cos(i*Pi)],
       r = 0 .. 1, 0.1), j=0..2, 1/5), i = 0 .. 1, 1/10)], symbol=diamond, shading=zhue);