Instead of the method used by dharr you can use Equate:
 

restart;
A:=Matrix([[2*x,z],[-2*x-2,x-2]]);
B:=Matrix([[4*y,2*y],[w,3]]);
Id:=Matrix([[1,0],[0,1]]);
Equate(A-B,Id);
solve(%);

For information on Equate do ?Equate.

See ? gridlines.
You can do:
 

plots:-display(p1,p2,gridlines);

In defining a procedure as you do for W, nothing is evaluated.
Then you give the parameters a, b, c, ... concrete values.
Then you define W once again in the exact same way (!), but as I said nothing will be evaluated.
You could do this using unapply the second time:
 

restart;
W := Un -> -1/4*(b*d - e*q)^2/(b^2*(e + b*Un^2)) - (2*b*c*m^2 + 1/4*q^2)*Un^2/b + c*Un^4
b := -0.0003070816340;
c := 0.4461893869;
d := 0.01142581922;
e := 0.0007675000601;
q := -0.003704049149;
m := 1.423206983;
eval(W); #The constants are NOT evaluated
D(W);    #The constants are NOT evaluated
W(x); #The constants ARE evaluated
W1:=unapply(W(x),x);
D(W1);

The package DirectSearch is available from the Maple Application Center:
https://www.maplesoft.com/Applications/

DirectSearch:-GlobalOptima(G(x,y),[x=0..1,y=0..1]); # Minimum
-DirectSearch:-GlobalOptima(-G(x,y),[x=0..1,y=0..1]); # Maximum

PDEtools:-dchange can be used:

restart;
e:=Int(1/(3*u^4 + u + 3),u);
res1:=PDEtools:-dchange({u=y*x^(1/3)},e,[y],params={x});
res2:=Int(x/(3*x^2*y^4+x*y+3*x^(2/3)),y); # Result from acer's use of IntegrationTools:-Change
simplify(diff(res1-res2,y)); #0

Here is an animation in gamma:

The code:
 

restart;
local gamma;
Digits:=15;
#gamma:=0.318;
alpha:=-1;
beta:=1;
delta:=0.1
omega:=1.4;
ics := {v(0) = 0, x(0) = 0};
sys := {diff(v(t), t) = -delta*v(t) - alpha*x(t) - beta*x(t)^3 + gamma*cos(omega*t), diff(x(t), t) = v(t)};
sol := dsolve(sys union ics, numeric, parameters=[gamma],maxfun = 0, abserr = 0.1e-17, relerr = 0.1e-12);
with(plots):
sol(parameters=[0.35]);
odeplot(sol,[t,x(t)],0..800,numpoints=10000,size=[1200,default]);
odeplot(sol,[x(t),v(t)],0..800,numpoints=10000,size=[1200,default]);
odeplot(sol,[x(t),v(t)],2200..3000,numpoints=10000,size=[1200,default]);
Q:=proc(gamma,{range::range:=2200..3000}) if not gamma::realcons then return 'procname'(_passed) end if;
  sol(parameters=[gamma]);
  plots:-odeplot(sol,[x(t),v(t)],range,numpoints=10000,size=[1200,default]);
end proc;
Q(0.318);
animate(Q,[gamma],gamma=0.318..0.35,frames=100,size=[1200,default]);

It makes some sense to me to plot [t,x(t).v(t)] too. This is just as converting the the existing nonautonomous system into an automous one in 3 variables.
The code for Q can just be revised in this way:
 

Q:=proc(gamma,{range::range:=2200..3000, scene::list:=[x(t),v(t)]}) 
  if not gamma::realcons then return 'procname'(_passed) end if;
  sol(parameters=[gamma]);
  plots:-odeplot(sol,scene,range,numpoints=10000,_rest);
end proc;
## and then animate in 3 dimensions:
animate(Q,[gamma,scene=[t,x(t),v(t)]],gamma=0.318..0.35,frames=100);

The animation:

Maybe a less confusing version: Narrow the range and make t the third variable instead of the first:
 

animate(Q,[gamma,scene=[x(t),v(t),t],range=2900..3000],gamma=0.318..0.35,frames=100);

When animating set the FPS (frames per second) to 1.

Using MultiSeries:-limit gives the same result for both versions:

 

restart;
res1:=limit(CylinderU(0,CylinderU(0,x)),x=0);
res2:=CylinderU(0,limit(CylinderU(0,x),x=0));
res1M:=MultiSeries:-limit(CylinderU(0,CylinderU(0,x)),x=0);
res2M:=CylinderU(0,MultiSeries:-limit(CylinderU(0,x),x=0));

res1<>res1M, but res2=res2M = res1M= CylinderU(0, 1/2*2^(3/4)*sqrt(Pi)/GAMMA(3/4)).

Further evidence can be obtained from the defining differential equation for CylinderU:
 

# y'' - (x^2/4 + a)*y = 0; a = 0 here
ode:=diff(y(x),x,x)-x^2/4*y(x)=0;
dsolve(ode); # Uses Bessel functions
y0:=CylinderU(0,0); 
y1:=eval(diff(CylinderU(0,x),x),x=0);
ics:=y(0)=y0,D(y)(0)=y1; # initial conditions
sol:=dsolve({ode,ics});
CU:=unapply(rhs(sol),x); # This is CylinderU(0,x)
limit(CU(CU(x)),x=0);
ev1:=evalf[15](%);
CU(limit(CU(x),x=0));
ev2:=evalf[15](%);
ev_res2:=evalf[15](res2);

In fact we can prove the identity of rhs(sol) and CylinderU(0,x):
 

B:=convert(CylinderU(0,x),Bessel); # Uses BesselJ only
B2:=convert(rhs(sol),BesselJ);
simplify(B-B2) ; # 0

There is really no need to simplify. dsolve converts floats ( i.e. numbers like 0.123 ) into rationals (in this case 123/1000).
That makes the result look complicated.
A simpler version of the solution (sol) is then provided by evalf(sol):
 

restart;
ode:=diff(f(x),x)=0.0768*f(x)^(2/3)-0.0102*f(x);# f(1)=59. 
sol:=dsolve({ode,f(1)=59});
evalf(sol); # Looks simpler
plot(rhs(sol),x=0..5000);

@dharr I think you are quite right:
 

restart;
ODE := diff(y(t), t) = -0.000032229*y(t)^2 + 0.065843*y(t) - 15.103;
RHS:=subs(y(t)=y,rhs(ODE));
z1,z2:=solve(RHS=0,y);
plot(RHS,y=z1-100..z2+100);
SOL:=dsolve({ODE,y(0)=z});
plot(eval(rhs(SOL),z=z1),t=0..800); # Constant (roughly)
plot(eval(rhs(SOL),z=z1+0.001),t=0..800);
plots:-animate(plot,[rhs(SOL),t=0..800],z=z1+0.001..z1+10,trace=5);

Enlarging the image provided I can see that the curve has to satisfy (t,y) = (0,449), thus this would be the curve:
 

plot(eval(rhs(SOL),z=449),t=0..350);

To get a direction field you can do
 

DEtools:-DEplot(ODE,y(t),t=0..350,[y(0)=449],y=0..2000);

The image is here:

Try this:
 

convert~(Res2,units,min);

Result:
Vector([65*Unit('min'), 70*Unit('min'), 75*Unit('min')])

restart;
pde:=diff(u(x, t), t, t) + 3 + 2*diff(u(x, t), t) + 4*t + x^2 + x^3/3 + diff(u(x, t), t, x, x) + diff(u(x, t), x, x, x, x) = x*t^2;
ode:=y(x)+diff(y(x),x$2)+x=sin(x);
##########################
p:=proc(eq::equation) local d,fu,res;
  if not has(eq,diff) then return eq end if;
  d:=indets(indets(eq,specfunc(diff)),function);
  fu:=op(remove(has,d,diff));
  res:=selectremove(has,(lhs-rhs)(eq),fu);
  res[1]=-res[2]  ## minus put on res[2] (see nm's correction below)
end proc:
  
   
p(pde);
p(ode);

You could do:

DateDifference(d1, d2, 'units' = mo);
round(%);

Have you tried option trace, as in p:=procedure(....) option trace; local ...; etc end proc.
As an example consider this recent example from MaplePrimes https://mapleprimes.com/questions/238264-Period-Unlimited-Decimal-Fraction

where  I have edited JAMET's procedure and here with option trace:
 

periode:=proc(r::rational) option trace; local a,b,c,f,i,k,l,p,q,s:  
  s:=convert(evalf(abs(r)-trunc(r),50),string):  
  if  s[1]="." then s:=s[2..length(s)] fi:  
  a:=numer(simplify(r)): 
  b:=denom(simplify(r)):  
  if  igcd(b,10)=1 then ## 1
    q:=0:       
    p:=1:      
    while (10^(p) mod b)<>1 do  
      p:=p+1 od:  
    else      
      f:=ifactors(b)[2]:
      k:=0:l:=0:
      for i to nops(f) do
        if f[i][1]=2 then k:=f[i][2] fi:
        if f[i][1]=5 then l:=f[i][2] fi: 
      od:
     c:=b/(2^k*5^l):
     q:=max(k,l):
     if c=1 then p:=0 else p:=op(2,periode(1/c)) fi:
   fi:
[q,p,[seq(parse(s[k]),k=1+q..q+p)]]
end:

It works Ok for some examples:
 

periode(2/3);
periode(50/7);

But not for all:
 

periode(1007/200035);

periode:=proc(r::rational) local a,b,c,f,i,k,l,p,q,s:  
  s:=convert(evalf(abs(r)-trunc(r),50),string):  
  if  s[1]="." then s:=s[2..length(s)] fi:  
  a:=numer(simplify(r)): 
  b:=denom(simplify(r)):  
  if  igcd(b,10)=1 then ## 1
    q:=0:       
    p:=1:      
    while (10^(p) mod b)<>1 do  
      p:=p+1 od:  
    else      
      f:=ifactors(b)[2]:
      k:=0:l:=0:
      for i to nops(f) do
        if f[i][1]=2 then k:=f[i][2] fi:
        if f[i][1]=5 then l:=f[i][2] fi: 
      od:
     c:=b/(2^k*5^l):
     q:=max(k,l):
     if c=1 then p:=0 else p:=op(2,periode(1/c)) fi:
   fi:
[q,p,[seq(parse(s[k]),k=1+q..q+p)]]
end:
periode(2/3);
periode(50/7);

Notice that I changed "couvert"  to "convert".

Note: This fails:

periode(1007/200035);

This appears to happen if q = 0.

To see that q = 0 insert a print statement before the parse statement.
 

s:=convert(evalf(abs(1007/200035)-trunc(1007/200035),50),string); 
q,p:=0, 1818;
[seq(parse(s[k]),k=1+q..q+p)];

Using allvalues on sol gives two results. odetest produces something not zero though.

restart;
sol:=y(x) = (exp(RootOf(-sin(x)*tanh(1/2*_Z+1/2*c__1)^2+sin(x)+exp(_Z)))+sin(x))/sin(x);
ode:=diff(y(x),x)-cot(x)*(y(x)^(1/2)-y(x)) = 0;

sol1,sol2:=allvalues(sol);
odetest(sol1,ode,y(x));
odetest(sol2,ode,y(x));