You didn't miss anything.  There are just two stationary points, of which one is a saddle point and the other is a local minimum.

> S := dsolve([Eq1, Eq2, Eq3a, R(0) = 49200, H(0) = 4879, 
    W(0) = 55800],[R(t),H(t),W(t)], numeric, 
    method = rosenbrock,output=Array([5*i $ i=0..8]));
  ExportMatrix("myfile.csv", S[2,1], target=delimited, 
    delimiter=",");

In the "while E > M do ...", E is supposed to be the absolute value of f at the last point
computed.  When that falls below M = 0.001, you stop your iteration.  Then you go and ask for a new A, B and C, but E is still the same (the value that made you stop before), so your while loop stops before it starts.  I think you'll want to put your

E := M + 1;

some place after the

while (A <> n) do

and before the

while E > M do

1) You have a parenthesis wrong in the formula: it should be

20 + 74 * exp(-ln(.4)*t/2)

a) Use solve.

b) Use plot.

c)Use limit.

Try labelledcontourplot in my Maple Advisor Database,

<http://www.math.ubc.ca/~israel/advisor>

It looks like fieldplot doesn't have an appropriate option, but you can duplicate the missing functionality  "by hand" using arrow.  Thus:

> myfieldplot:= proc(f::list, r1::name=range, r2::name=range)
    local x,y,a,b,c,d,g,dx,dy;
    if not typematch(r1, x::name=a::realcons .. b::realcons) then
      error("expected name = (real constant) .. (real constant), got  ",r1)
    end if;
    if not typematch(r2, y::name =c::realcons .. d::realcons) then
      error("expected name = (real constant) .. (real constant), got  ",r1)
    end if;
    if not hasoption([args[4..-1]],grid=[posint,posint],g) 
       then g:= [10,10] 
    end if;
    dx:= (b-a)/g[1]; dy:= (d-c)/g[2];
    plots[arrow]([seq(seq(eval([[x,y],f],{x=a+i*dx,y=c+j*dy}),
       i=0..g[1]),j=0..g[2])],shape=arrow);
 end proc;

This was written to accept only the grid option, but others can be included if desired.

> plots[animate](myfieldplot,[[t/10,t/10], x=-10..10,y=-10..10], 
      t = -10 .. 10);

Assign it to a variable first.

> P := plot(...):
  P;

[ shows the plot on the screen ]

> plotsetup(...);
  P;
  plotsetup(default);

[ writes it to a file]

The best I can suggest for automatic simplification is

> simplify(%, size);

 ((-f(zn)+D(f)(zn)*zn)*f(xn)^2+(f(zn)^2-2*f(zn)*D(f)(zn)*zn)*f(xn)+D(f)(zn)*xn*f(zn)^2)/D(f)(zn)/(f(xn)-f(zn))^2

It's better to use Matrix rather than the old, deprecated matrix.

> F:= <<0 | 1>, <0 | 0>>;
  C:= <0, c>;
  G:= <1 | 0>;
  d := m;

(I'm using d instead of D, because D is reserved for the differentiation operator).  Note that d is a scalar rather than a Vector or Matrix.

> S := Matrix(2,2,(i,j) -> s[i,j](t));
  R:= map(diff,S,t) - (F.S + S.F^%T - S.G^%T * (d . d)^(-1) . G . S + C . C^%T);
  DEs := {seq(seq(R[i,j], i=1..2), j=1..2)};
   

Unfortunately the symbolic dsolve can't do much with this system, so we'll try numeric.  Your initial conditions would produce a trivial solution, so I'll make it slightly more interesting.

> Res:=dsolve(subs(c=0,m=1,DEs) union {seq(seq(s[i,j](0)=i,i=1..2),j=1..2)}, 
     numeric, output=procedurelist);
   Sol:= t -> subs(Res(t), S);
   Sol(5); 

[ .520368044229340798    .0867280208733355874 ]
[ .0946123788437747138   .0157687304230208666 ]

> plots[odeplot](Res, [seq(seq([t,s[i,j](t)],i=1..2),j=1..2)],0..5, 
     legend=[seq(seq(s[i,j](t),i=1..2),j=1..2)]);

You could use the circle command in the plottools package, together with display from the plots package.  See the help page

?plottools,circle

Equivalent is a command in the Logic package.

In that package, the logical operators are &and, &or, &not, &implies etc.

So for example:

> with(Logic):
  Equivalent((p &implies q) &and (q &implies r), 
        p &implies r);

   false

(no, those expressions are not equivalent)

> Equivalent((p &implies q) &and (q &implies r), 
       ((p &or q) &implies (q &and r)));

  true

(yes, those are equivalent).

By the way, it is easy to post your Maple worksheet.  Click the green up-arrow button and follow the directions.

AFAIK there isn't anything specific for delay differential equations, but you can just integrate interval-by-interval.  Let me do a slightly simpler example:
 

y'(t) = 1 + y(t-1) - y(t-2), with y(t) = t^2 for t <= 0. 

> y[-1]:= t -> t^2;  # this is y(t) for t <= -1
    y[0]:= t -> t^2;    # this is y(t) for t <= 0 
    for j from 1 to 10 do
        y[j]:= unapply(y[j-1](j-1) + int(1 + y[j-1](s-1) - y[j-2](s-2), s = j-1 .. t), t); # this is y(t) for j-1 <= t <= j  
   end do;
    plot(t -> y[ceil(t)](t), -2 .. 10);

> with(LinearAlgebra):
   A:= Matrix([[1,1,1,1,1,1],[1,0,-1,1,0,-1],[2,1,0,2,1,0],[4,1,2,0,4,1],[1,-2,1,-6,0,1],[5,4,-1,-2,5,0]]);
   b:= <1,0,1,0,1,0>;
   LinearSolve(A,b);

< 3/20*_t[3]+11/12, 13/20*_t[3]+5/12, _t[3], -1/5*_t[3], -11/20*_t[3]-5/4, -21/20*_t[3]+11/12 >

 There are solutions, so it is consistent.

The most usual sufficient condition is Fubini's theorem, of which one case is: if f(x,y) is a measurable function on (a,b) x (c,d) (with some of a,b,c,d allowed to be infinite) and int(int(abs(f(x,y)), x=a..b),y=c..d) < infinity, then int(int(f(x,y), x=a..b), y=c .. d) = int(int(f(x,y), y=c..d), x=a..b).  In your case this will do it for you as long as p > 0 and q is real, because |exp(-p*x) * cos(q*x)| <= exp(-p*x), and int(exp(-p*x), x=0..infinity) < infinity.

The Hessian determinant is not enough.  You need to check whether the Hessian matrix is positive semidefinite.  There are better ways of doing that than computing eigenvalues, though.  In Maple, you can use the IsDefinite function in the LinearAlgebra package.

Determining whether a function is always nonnegative on a given domain is not an easy task in general, but "is" can sometimes do it.

For example:

> F:=  x^3 + x*y^3 - x*y;
  Q:= LinearAlgebra[IsDefinite](VectorCalculus[Hessian](F, [x,y]), query=positive_semidefinite);

Q := 0 <= -9*y^4+6*y^2-1+36*x^2*y and 0 <= 6*x*y+6*x

> is(Q) assuming  x > 1, y > 1/2, y < 1;

       true