They are mathematically equivalent, but `simplify/D` does not implement this. But,

e1 := (D[1]@D[1])(f)(x);
convert(e1,diff): convert(%,D): lprint(%);
    ((D@@2)(f))(x)

p1:=plot3d([[x,y,-y^2],[x,y,x^2]], x=0..1,y=0..x):
p2:=plot3d([[1,y,z]], z=-y^2..1,y=0..1,transparency=0.5):
p3:=plot3d([[x,0,z]], x=0..1,z=0..x^2,transparency=0.5):
p4:=plot3d([[x,x,z]], x=0..1,z=-x^2..x^2,transparency=0.5):
plots[display](p1,p2,p3,p4,labels=["x","y","z"]);

# less memory used

myPDF := (v,t) -> piecewise(t < 0, 0, t < v, 1/v, 0):
myD:= Distribution(PDF = curry(myPDF, 1)):  #v=1
Mean(myD);
    1/2

https://oeis.org/A007018

J:=Int( exp(-I*k*x)/cosh(x), x=-infinity..infinity):
value(IntegrationTools:-Change(J,exp(x)=t));

1. Due to the wonderful 2D math, a variable named _ appeared. This must be corrected.

2. The equation for solve is too complicated to be solved symbolically.

3. For fsolve you must have numerical values for parameters.

PW1:=proc(e::specfunc(anything,piecewise))
local f, o:=op(e);
f:=proc()
  if   nargs=3 then '`if`'(args)
  elif nargs=2 then '`if`'(args,0)
  else '`if`'(args[1],args[2],f(args[3..-1])) fi
end;
f(o)
end:

PW:=e->subsindets(e, specfunc(anything, piecewise), PW1):
####################

PW(piecewise(x<1, 10, x<2, 20, x<3, 30));
    `if`(x < 1, 10, `if`(x < 2, 20, `if`(x < 3, 30, 0)))

f:=piecewise(x < 0, piecewise(y < 0, 0, 1), piecewise(y < 0, 1, 0)):
PW(f);
    `if`(x < 0, `if`(y < 0, 0, 1), `if`(y < 0, 1, 0))

L:=(u,v)->u(v(f_(x,t)))-v(u(f_(x,t))):
alias(``=f_(x,t));
#########################
v:=f->diff(f,x):
w:=f->x*diff(f,t):
L(v,w);
   

u:=f->x*t*diff(f,x,t):
L(u,w):  simplify(%);

   

I was curious to see how much the compiler can help in such problems.

It seems that in Carl's version (which is very fast) the compiler is not very useful.
I used practically my original version, made compilable.

Here it is:

Tri1:=proc(n::integer,V::Vector(datatype= float[8]),
     X::Vector(datatype= float[8]),Y::Vector(datatype= float[8]),
     Ax::float[8],Ay::float[8],Bx::float[8],By::float[8],Cx::float[8],Cy::float[8])
   local k::integer,a::float[8],b::float[8];
   for k to n do
     a:=V[k]; b:=V[n+k];
     if a+b>1 then a:=1-a;b:=1-b fi;
     X[k]:=Ax+(Bx-Ax)*a+(Cx-Ax)*b;
     Y[k]:=Ay+(By-Ay)*a+(Cy-Ay)*b
   od;  
end:

CTri1:=Compiler:-Compile(Tri1):

Tri:=proc(n,Ax,Ay,Bx,By,Cx,Cy)
local
  V:=LinearAlgebra:-RandomVector(2*n, generator= 0..1., datatype= float[8]),
  X:=Vector(n,datatype= float[8]), Y:=Vector(n,datatype= float[8]);
global Ctri1;
  CTri1(n,V,X,Y,Ax,Ay,Bx,By,Cx,Cy);
X,Y
end:
####
Ax:=-1: Ay:=0:  # Triangle T = ABC
Bx:=0:  By:=4:
Cx:=2:  Cy:=1:
n:=6000:           # number of random points
T:=CodeTools:-Usage(Tri(n,Ax,Ay,Bx,By,Cx,Cy),iterations=100):
memory used=195.55KiB, alloc change=18.38MiB, cpu time=940.00us, real time=920.00us, gc time=0ns

Compare with Carl's version (not compiled)

T:= CodeTools:-Usage(SampleTriangle([[-1,0],[0,4],[2,1]], 6000), iterations=100):
memory used=0.86MiB, alloc change=8.21MiB, cpu time=24.18ms, real time=6.52ms, gc time=1.25ms

Note that I have used iterations=100, otherwise the timing was not stable enough.

The definition is simple and clearly presented in the help page.

HeunT is the solution of the differential equation (depending on three parameters):

f(0)=1,  f'(0)=0

It is an entire function. For other properties:

FunctionAdvisor(HeunT);

Yes, it seems that the numerical algorithm in sum fails here.

But we can verify manually (after a simple estimation of the remainder):

evalf[600]( add( floor((exp(Pi)-Pi)*n)/3^n, n=1..5000) ):
evalf(evalf[600](%-29/2));
     -4.142437940*10^(-531)

Ax:=-1: Ay:=0:  # Triangle T = ABC
Bx:=0:  By:=4:
Cx:=2:  Cy:=1:
n:=4000:        # number of random points

a:=Vector(n): b:=Vector(n):
r:=rand(0.0 .. 1.0):
for i to n do
  x:=r(): y:=r():
  if x+y>1 then  x,y := 1-x,1-y  fi;
  a[i]:=x: b[i]:=y:
od:
u:=Vector(n,1):

plot(
Ax*u+(Bx-Ax)*a+(Cx-Ax)*b, Ay*u+(By-Ay)*a+(Cy-Ay)*b,
symbol= point,style=point,scaling=constrained);

# The distribution is uniform i.e.
# Prob(X in D) = Area(D)/Area(T).

n:=8000:   #ellipse - nonuniform distribution (naive attempt)
a:=Vector(n,i->r()):
b:=Vector(n,i->r()):
plot(3*a*~cos~(2*Pi*b),2*a*~sin~(2*Pi*b),
symbol= point,style=point,scaling=constrained);

a:=sqrt~(a):   #ellipse - uniform distribution
plot(3*a*~cos~(2*Pi*b),2*a*~sin~(2*Pi*b),
symbol= point,style=point,scaling=constrained);

If you can't stop the executecution of your worksheet, use mine:

x:=0.0:
for i to 10^8 do x:=x+1/i od:x;
###
for i to 10^8 do x:=x+1/i od:x;
###
for i to 10^8 do x:=x+1/i od:x;

which can be interrupted. :-)

The wmf format supports both vector and bitmap graphics, but unfortunately Maple produces only bitmap.

For Windows 64bit (at least), there is no support for vector 3D graphics! (postscript is supported but totally wrong!). I use Maxima for this purpose.

extrema(y,{},{x,n});

or

yy:=subs(x=t/n,y);
extrema(yy,{},t);