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