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));

I did it. In order to use it you have to agree to the terms of use.
Wanting to test it I agreed to the terms.  After the testing I took that back. That is possible in the same place.

If you want to avoid the weird looks of the definition ot totvol in your second run, use unapply instead of what you got:
 

totvol := unapply(Pi*0.4^2*ht^2*ht + 1/3*Pi*0.4^2*ht^2*2/3*ht,ht);

The output is simple:  totvol := ht -> 0.1955555556*Pi*ht^3

Illustrating vv's point I tried starting solutions at x=1 for several values of y(1) and plotting the results:
 

restart;
de1 := diff(y(x), x) = y(x)/x - 5^(-y(x)/x);
S:=seq(rhs(dsolve({de1,y(1)=k})),k=0..20);
plot([S],x=0..1);
plots:-display(seq(plot(S[k],x=0..1),k=1..21),insequence);# An animation

The animation is attached: