dN/dt is the rate of change of the number of people infected.  This is supposed to be proportional to N*(P-N).  Do you know what proportional means?

Note that like all inverse trig functions, arctanh is a multi-valued function.  Its branch points are at 1 and -1.   Since tanh(w + I*Pi) = tanh(w), if w is a possible value of arctanh(z) then so is w + I*Pi.  Thus if there is any solution of tanh(w) = z, there is one with imaginary part in the interval [-Pi/2, Pi/2].  Since tanh(t + I*Pi/2) = tanh(t - I*Pi/2) = coth(t) goes from +infinity to +1 as t goes from 0 to infinity and from -infinity to -1 as t goes from 0 to -infinity,
for every real z with |z| > 1 there is a value of arctan(z) with imaginary part +Pi/2 and another with imaginary part -Pi/2.  The convention Maple uses is that the imaginary part is
-Pi/2 for z > 1 and +Pi/2 for z < -1, and in the interval (-Pi/2, Pi/2)  for all other z. 

You seem to have left something out of your definition of p.  Also, is pn the same as p?

The most challenging part of this, I think, is to have rotated characters.  As far as I know there is no option for this: characters are always plotted in their standard orientation (except that axis labels can be made to be vertical rather than horizontal using the labeldirections option).  So you'd have to actually draw the characters (using vector-graphics definitions of them) rather than using something like textplot.

If the population is P, and N people are infected, then what is the product of the number infected and the number not infected?

Yes, it's certainly a bug.

You probably intended to have r > 0.  With that assumption, it works.

> assume(r>0);
    simplify(asympt(B(n,1)/B(n,0)-1,n,3));

Compare (with this assumption on r)

> asympt(r^(n-1), n);

(correct)

> asympt(q^(n-1),n); # no assumption on q

(wrong)

No, pdsolve/numeric is still only for "time-based" pde's.

b:= matrixplot(a,heights=histogram,colour=((x,y) -> a[floor(x),floor(y)])):
plottools[transform]((x,y,z) -> [(50+z)*cos(x*Pi/6),(50+z)*sin(x*Pi/6),y])(b);

Do you mean something like this?

> plottools[transform]((x,y,z) -> [(50+z)*cos(x*Pi/6),(50+z)*sin(x*Pi/6),y])(b);

surfdata is probably what you want if the curves have the same number of points.  Otherwise you might look at this discussion from last year.

> PDE:= 2*diff(U(x,y),x,x) + 6*diff(U(x,y),x,y) + 3*diff(U(x,y),x) + 2*U(x,y) = 0;
   seq(frontend(coeff,[PDE,t]),t=[diff(U(x,y),x,x),diff(U(x,y),x,y),diff(U(x,y),y,y),diff(U(x,y),x),diff(U(x,y),y),U(x,y)]);

                           2, 6, 0, 3, 2, 0

Are those "-" actually in the file, or did you just put them there for spacing?
Are the numbers in the file all positive integers?  What about values other than 2 and 768?  Is there a way to translate from your four numbers to a direction vector?

It's really hard to say without seeing the actual worksheet.  You can upload it using the green up-arrow button.

One possible cause of this sort of behaviour is that something connected to your equation is assigned a value in the worksheet where it doesn't work.   We just had a case similar to that in this thread.

Are you talking about the sphere command in the geom3d package, or the one in the plottools package?  One of these might be useful for (a), the other for (b).

The problem seems to be the fact that x[0] has been assigned a value.  It often causes trouble if you try to use a name (x in this case) as a symbolic variable when an indexed version of that variable has been assigned a value (which means your variable is a table rather than an unassigned name).  Change x[0] to x0 or X[0] and it should all work fine.