Now the dictionary appears as "Definition".

Looking for ?Definition  the first definitions appear in the search list, but not all, probably because the length of the list is limited.
But searching e.g. ?Definition,taylor,  the  definitions related to taylor appear.
In my opinion, the organization and the seach system used to be better in earlier versions.

I think that a good idea is use tables with multiple indices.
E.g. instead of
A :=  table([x=table([u=1]), y=table([v=table([p=2]), w=3])]);
use:

A[_root]:=[x,y];
A[x]:=[u];
A[x,u]:=1;
A[y]:=[v,w];
A[y,v]:=[p];
A[y,v,p]:=2;
A[y,w]:=3;

So, the leaves will be:

select( t -> not(type(t,list)), [entries(A,nolist)]);
And the corresponding path is in
indices(A);

 

Your equation is:

tan(Pi*t) = tan(2*Pi*t^2);
Hence, the solutions  are

1/4+(1/4)*sqrt(1 - 8*_Z), 1/4-(1/4)*sqrt(1 - 8*_Z);

(_Z integer; <=0 for real solutions).

The best way would be to save the generated procedures.

restart;
n := 40: h := 40:
a := h*i/(2*n):

sys1 := [diff(x(t), t) = 2*t+1, diff(y(t), t) = 1+5*t, x(0) = 1, y(1) = 0]:
for i from 0 by .5 to 1 do res[i] := dsolve(sys1, numeric) end do:
save res, "d:/mymaplefiles/odes.txt":
# ...
restart;
read "d:/mymaplefiles/odes.txt":
res[0](1/2);

Actually in 1D math one may use any number n of dots (n>1) instead of n=2.
So,  seq(2..........7)  is OK.
The situation in 2D math is the same, unless there is no space after the last dot . So, seq(2.......  7) is OK. Without spaces, the number of dots must be even (otherwise the last dot is considered as a decimal point). It is interesting.
 

You are missing a multiplication sign in z[k+1]  after (1 + 1/beta).

2D input math in action again!

Generate each plot in a different color.
For example,

col := k -> COLOR(HUE, k/(2*nops(n)+1)):

#  inside the loop:
# for ...
fplt[j] := plots[odeplot](sol1, [eta, diff(diff(f(eta), eta), eta)], axes = boxed, color=col(j));
tplt[j] := plots[odeplot](sol1, [[eta, theta(eta)]], axes = boxed, color=col(j+n))
# ...end do;

In Windows insert in your code something like

system("C:/Program Files/Windows Media Player/wmplayer.exe  d:/yourdir/yourmelody.mp3");

I think I found the location where the random behavior resides.
It is the indices(T) where T is a table containing local variables in indices.
That's all we can do because indices is builtin.

p:=proc()
local lambda,T;
T[lambda-1]:=1; T[lambda-2]:=1; T[lambda-3]:=1;
indices(T);
end:

p();


At each execution, the order changes.
Note that for a nonlocal lambda, everything is OK.

foldr and foldl are simply executing a loop. Why not doing it yourself:

p:=1:
for x in X do p:=mods(mods(x,m)*p, m) od;

Maybe mods(x,m)  could be eliminated and replaced by x;  this depends on the magnitude of your numbers.

For sin(t)<0  e.g. t in =-Pi .. 0 you have a simple ODE x' = - B*x
for which an initial condition is not provided.

For x in 0 .. Pi it results an integro-differential equation in x(t), but probably Maple can't solve it (numerically or symbolically) and some special approximation method will be needed.

This  solution could maybe used as a starting approx.

restart;
Digits:=100;
EQ1 := -3.000000000*10^(-8)+3.815358072*sin(3.141592654*10^7*x)+9.534375000*10^(-30)*(diff(w(x), x, x, x, x))-2.383593750*10^(-60)*(diff(w(x), x, x, x, x, x, x))-5.085000000*10^(-13)*(diff(w(x), x))*(diff(u(x), x, x))-7.627500000*10^(-13)*(diff(w(x), x))^2*(diff(w(x), x, x))-5.085000000*10^(-13)*(diff(w(x), x, x))*(diff(u(x), x))+0.2410290000e-5*(diff(w(x), x, x));
EQ2 := 5.650000000*10^(-20)*(diff(u(x), x, x, x, x))-226000000000*(diff(u(x), x, x))-226000000000*(diff(w(x), x))*(diff(w(x), x, x));
bc:= {u(0) = 0, u(L) = 0, w(0) = 0, w(L) = 0, (D(u))(0) = 0, (D(u))(L) = 0, ((D@@2)(w))(0) = 0, ((D@@2)(w))(L) = 0, ((D@@4)(w))(0) = 0, ((D@@4)(w))(L) = 0};
Order:=8;
s:=dsolve( {EQ1,EQ2}, {u(x),w(x)}, 'type=series'):
s:=evalf(convert(s,polynom));
s1:=eval(s,bc);
#indets(%,function);
substs:=[(D(w))(0)=w0, ((D@@2)(u))(0)=u2, ((D@@3)(u))(0)=u3, ((D@@3)(w))(0)=w3, ((D@@5)(w))(0)=w5];
s2:=eval(s1,substs);
uf:=unapply( eval(u(x),s2), x);
wf:=unapply( eval(w(x),s2), x);
subs([u=uf,w=wf],bc);
sys:=eval(%) minus {0.=0};
#nops(sys);
L := 100.*10^(-9);
ok:=solve(sys);
solux:=eval(uf(x),ok[1]);
solwx:=eval(wf(x),ok[1]);
plot(solux,x=0..1);
plot(solwx,x=0..1);

The graph of u:

For x = 0 .. L it is:

You cannot, Maple does not know about it.

But usually it is straightforward (except pathological examples) to express a Stieltjes integral as a Riemann one.
Note that Maple also knows improper Riemann integrals such as int(1/sqrt, 0..1).

A recursive version, simple and fast.
It returns a sequence of lists.

In my tests it was the fastest for (n,k) = (22,11).

rcomb:=proc(n,k,p:=[])
  if k=n then return [seq(1..n),op(p)] fi;
  if k=1 then return seq([i,op(p)],i=1..n) fi;
  rcomb(n-1,k,p),rcomb(n-1,k-1,[n,op(p)])
end;

 

nops([rcomb(20,5)])=binomial(20,5);
                         15504 = 15504

rcomb(5,3);
    [1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4], [1, 2, 5], [1, 3, 5], [2, 3, 5], [1, 4, 5], [2, 4, 5], [3, 4, 5]

Edit. Corrected code. Unfortunately some speed was lost, but still decent.

The problem is nonlinear and nonconvex-nonconcave, so global min/max are not easy.
Your V must be a row rector, but I will use columns.


restart;
with(LinearAlgebra):
n:=2;
A:=RandomMatrix(n);
B:=RandomMatrix(n);
V:=Vector(n,symbol=v);
f := (V^+ . A^+ . A . V) / (V^+ . B^+ . B . V);
# dom := seq( v[i]=-1..1, i=1..n),location;
dom := location:
maximize(f, dom);
evalf(%);
minimize(f, dom);
evalf(%);

Note that for n>3 this will not work because the exact solutions will need the eigenvalues of a matrix of order n.

Mathematically this problem should be approached differently:
the objective function is the square of  Norm(A.V)/Norm(B.V)  and this
reduces to compute the norm of a matrix  T^(-1).A.T  where T is an orthogonal matrix deduced from B.

Edit. Minimize works but gives only approximate (float) results and the min/max could be only local.

n:=6;
dom1 := seq( v[i]<=1, i=1..n),  seq( v[i]>=-1, i=1..n):
A,B:='LinearAlgebra:-RandomMatrix(n)'$2;
V:=Vector(n,symbol=v):
f := (V^+ . A^+ . A . V) / (V^+ . B^+ . B . V):
Optimization[Minimize](f, {dom1});