There is no such thing as ODEplot.  There is odeplot in the plots package, and there is DEplot in the DEtools package.  I suspect the latter may be more suitable for your needs.

> with(DEtools): with(plots):
> de:= D(y)(x) = (x* (1 - .5*x - .5*y(x)) ) / (y(x)* (-.25 + .5*x));
  DEplot(de,y(x),x=0.5 ..3,{[y(1)=2],[y(1)=4]},arrows=none,linecolour=[red,blue]);

Now I don't know what you mean by the "critical points" of the solutions in this case.  I suspect you don't really mean the solutions of this differential equation, which don't really have critical points, but rather trajectories of the related system
 

> sys:= {D(y)(t) = x(t) * (1 - .5*x(t) - .5*y(t)), D(x)(t) = y(t) * (-.25 + .5*x(t))};

Perhaps something like this?

> DEplot(sys,[x(t),y(t)],t=-5 .. 10, [[x(0)=1,y(0)=2],[x(0)=1,y(0)=4]],arrows=none,
linecolour=[red,blue],stepsize=0.02); P1:= %:
  sols:= [seq(dsolve(sys union {x(0)=1,y(0)=yn},numeric,output=listprocedure),yn=[2,4])];
  ts:= [seq(fsolve(subs(sols[i],y(t)),0..5), i=1..2)];
  pts:= [seq([subs(sols[i],x(t))(ts[i]),0],i=1..2)];
  P2:= pointplot(pts,symbol=solidcircle,symbolsize=16,colour=black):
  P3:= arrow([pts[1]+[1.1,1.1],pts[2]+[1.1,1.1]],[<-1,-1>,<-1,-1>]):
  display([P1,P2,P3],scaling=constrained,view=[0..4,-1..4]);

For some strange reason, the Input format sometimes changes from Filtered HTML (which is what it should be) to Plain Text.  Click on "Input format" and you can change it. 

Actually sets would not be appropriate, because the numbers of repeated elements must be taken into account.  Lists would be better.  But for the sake of efficiency, a 1D Array might be better still.

> use ArrayTools in 
     IsEqual(sort(Reshape(A,[n^2])), sort(Reshape(B,[n^2])))
   end use;
   

The semicolon/colon that affects printing here is the one after the end if.  When you have a statement (such as an if statement)
that contains statements within it, what determines whether output is printed is the semicolon or colon at the end of the overall (top-level) statement.  If you want the graphs shown, use print.  Thus

if
...
 print(display(gArr));
...
end if:

1) You could use the geometry package, but an algebraic approach might be simpler.

2)  It's really hard to give a hint on something so easy.  Are you allowed to use convert?

Actually |z+2i|=2 would correspond to x^2 + (y+2)^2 = 4, not 2.

> with(plots):
   Aplot:= implicitplot(max(x^2+(y+2)^2-4,1-x^2 - (y+2)^2) <= 0, x = -2 ..2, y = -4 .. 0,   
      grid=[50,50],filledregions=true,gridrefine=3,crossingrefine=3,transparency=1/2):
   Bplot:= implicitplot(max(abs(argument(x+(2+y)*I))-Pi/3,1-abs(argument(x+(2+y)*I))),
     x=-2..2,y=-4..0, grid=[100,100],filledregions=true,gridrefine=3,crossingrefine=3,
     transparency=1/2, coloring=[green,white]):
  display(Aplot,Bplot);

I would interpret "compute the maxiumum numerical number" as: take the maximum of those members of the list that are of type numeric.
Thus:

> MY_MAX:= proc(L::list)
  local Ln;
  Ln:= select(type,L,numeric);
  if Ln =[] then ERROR
  else max(Ln)
  end if
 end proc;

And for the other question:

> PRIME_BRACKET :=  nextprime @ (t -> t-1) @ ceil;

pointplot in the plots package accepts a more efficient data structure, an Array of type float[8], but it "unpacks" that into an expression sequence of lists.  You might use quickplot in my Maple Advisor Database, <http://www.math.ubc.ca/~israel/advisor>.  Once that is installed, you could
try something like this:

> a:= 1.2; b:= 0.4; imax:= 100000;
   XY:= [seq(Matrix(imax+1,2,datatype=float[8],order=C_order),i=1..20)]:
   P0:= <0.1|0>;
   for N from 1 to 20 do
     XY[N][1,1..2]:= P0;
     for i from 1 to imax do
       XY[N][i+1,1..2]:= <1+XY[N][i,2] - a*XY[N][i,1]^2 | b*XY[N][i,1]>
     od;
     P0:= XY[N][imax+1,1..2];
     P[N]:= quickplot(XY[N],style=point,symbol=solidcircle,symbolsize=4,colour=black,
       axes=boxed,font=[arial,normal,18]):
   od:

The result doesn't seem too bad: each P[i] has length only 400146.  The Maple 12 Java interface crashed for me with an animation produced using all these plots, but Maple 13 (on Windows Vista) manages to do it.

> plots[display]([seq(P[i],i=1..20)],insequence=true);

Hint: try to keep things simple.  How about a diagonal matrix?

I suspect you're using 2D math input and have a space between the ln and the (n).  That causes implied multiplication to kick in, so Maple interprets the input as

> q := seq([ln*(n)/n, a(n)], n = 1 .. 25);

As the error message says: dsolve expects to get an equation, or a set or list of equations.  What you gave it was a list of two lists
(since sys and ics are both lists).  You could try

> ivpsolution := dsolve([op(sys),op(ics)]);  

However, the result (if Maple can find it at all)  would be very unpleasant: this is trying to solve a 5 x 5 linear system symbolically, and the solution will involve the exact roots of a quartic polynomial.  What happened when I tried it (in Maple 12 Classic under Windows Vista) was that it eventually crashed Maple.  You might prefer doing it with numerical values for those roots, something like this:

> vars:= [AT(t), TB(t), TD(t), MB(t), MD(t)];
  A:= LinearAlgebra[GenerateMatrix](map(rhs,sys), vars)[1];
  LinearAlgebra:-MatrixExponential(evalf(A),t) . Vector(map(rhs,ics));
  ivpsolution:= zip(`=`,vars,convert(%,list));

ivpsolution := [AT(t) = 591.7072634*exp(-.4894461536e-11*t)+128.7002551*exp(-.8678791583e-1*t)+117.7938254*exp(-.1522752038*t)-99.44158106*exp(-5.010041092*t)-48.7164273*exp(-.2356659624e-1*t), TB(t) = 281.3703586*exp(-.4894461536e-11*t)+328.3378355*exp(-.8678791583e-1*t)-131.0058251*exp(-.1522752038*t)+1.028707853*exp(-5.010041092*t)-29.73020094*exp(-.2356659624e-1*t), TD(t) = 437.6862264*exp(-.4894461536e-11*t)-357.1502637*exp(-.8678791583e-1*t)+62.4403426*exp(-.1522752038*t)-.1148907778e-1*exp(-5.010041092*t)-135.9648030*exp(-.2356659624e-1*t), MB(t) = 8.213515288*exp(-.4894461536e-11*t)+1.805205553*exp(-.8678791583e-1*t)+1.695193504*exp(-.1522752038*t)+688.9365716*exp(-5.010041092*t)-.650168052*exp(-.2356659624e-1*t), MD(t) = 3528.045758*exp(-.4894461536e-11*t)-101.6561011*exp(-.8678791583e-1*t)-50.96640480*exp(-.1522752038*t)-590.5122088*exp(-5.010041092*t)+214.9245629*exp(-.2356659624e-1*t)]
  

Note that the -.4894461536e-11*t should actually be 0: obviously there has been some roundoff error.  To fix that:

> ivpsolution:= fnormal(ivpsolution);

Well, you could do something like this.  First give values to p and lambda.  Then start with x[1] = 1 with probability p or 0 with probability 1-p, and get each x[j] from the previous one using the relevant conditional probabilities.

> p:= 0.3; lambda:= 0.4;
   P((X[i]=1) &given (X[i-1]=0)) := (1-lambda)*p/(1-p);
   P((X[i]=1) &given (X[i-1]=1)) := lambda;

> with(Statistics):
   U:= RandomVariable(Uniform(0,1)):
   S:= Sample(U,1000):
   if S[1] <= p then x[1]:= 1 else x[1]:= 0 end if:
   for j from 2 to 1000 do
      if S[j] <= P((X[i]=1) &given (X[i-1]=x[j-1])) then x[j]:= 1 else x[j]:= 0 end if
   end do:

One thing I might note is that your f is an expression, not a function, so IBC[i] should use p(a,i*T)=f, not f(a).  However, that isn't the source of the problem.  Your PDE is diff(p(a, t), a)+diff(p(a, t), t) = -p(a, t) which is stable in the forward direction but not in the backward direction.  pds[i+1] has its initial condition at time t=i*T, but you're trying to plot the solution from t=T to t=(i+1)*T.  If you have q[i+1] use t = i*T .. (i+1)*T, it should be OK.

> P(X[i]=1):= p; # assumed
   P(X[i] = 0):= 1 - P(X[i]=1);
   P((X[i] = 1) &given (X[i-1] = 1)) := lambda; # assumed
   P((X[i] = 0) &given (X[i-1] = 1)):= 1 - P((X[i] = 1) &given (X[i-1] = 1));
   

P(`&given`(X[i] = 0, X[i-1] = 1)) := 1-lambda

> solve({seq(P(X[i] = j) = add(P((X[i] = j) &given (X[i-1]=k)) * P(X[i]=k), k=0..1), j=0..1)},
  {seq(P((X[i]= j) &given (X[i-1]=0)), j=0..1)});

{P(`&given`(X[i] = 0, X[i-1] = 0)) = -(1-2*p+lambda*p)/(-1+p), P(`&given`(X[i] = 1, X[i-1] = 0)) = (-1+lambda)*p/(-1+p)}

What you want is a parametric plot: something like

> plot([ g, x^2, x = -5 .. 5]);

See the help page ?plot,parametric