Try this:

restart;
eta:= (x,y)-> y*sqrt(U*v/x);
alias(eta= eta(x,y));
Phi:= (x,y)-> sqrt(U*v*x)*f(eta);
e3:= vx = diff(Phi(x,y),x);
e4:= vy = diff(Phi(x,y),y);

sort([B||(1..3)], (A,B)-> (A[2]-A[1]) > (B[2]-B[1]))[1];

This looks at the width of the interval, not just the upper limit of the interval.

The sort is builtin and very fast. But if you were doing a very large number of intervals (millions probably), it may be faster to code a linear search for the max.

It's very easy:

L:= [[0,2,3,0,0], [2,3,0,0,0], [1,1,2,2,0]]:
add(1/mul(x!, x= LL), LL= L);

Also, it is not necessary that the lists be the same length. So there is no need to pad with zeros.

Clearly there is a bug here. Here is a workaround: Bring the w inside and switch the order of summation:

sum(sum(binomial(2*q,n-1)*z^q*w^n/q!, n=1..infinity),q=0..infinity);

Thanks for a quite interesting problem.

Unfortunately, dsolve doesn't allow explicit parameters for BVPs. But we can trick it into taking them by calling dsolve within a procedure whose parameters are the parameters that we want. I take the first three of your four bc and use them in dsolve. Then I make an implicitplot using the fourth bc.

Eq:= diff(f(y),y$3) + f(y)*diff(f(y),y$2) - diff(f(y),y$1)^2 -
     S*(diff(f(y),y$1) + y/2*diff(f(y),y$2))=0:
bc:= D(f)(0)=1, f(0)=0, (D@@2)(f)(b)=0:

F:= proc(_b, _S)
local Sol:= dsolve(eval({Eq,bc}, {S= _S, b= _b}), f(y), numeric);
     eval(f(y), Sol(_b)) - _S*_b/2
end proc:
     
plots:-implicitplot(F, 1..3, 1..3);

It will take some guess work adjusting the ranges to get the plot that you want. There are numerous options to implicitplot that can be used for refinement. A plot3d of F may help you find the level curve F=0. Note that b is the first argument of F, so b is on the horizontal axis.

My oh my, you've certainly made this far more complicated and redundant than it needs to be. You don't need odeplot, pointplot, or display; simply plot will do. Simplify your code to this:

restart:
A := 0.2e-1: B := 10^(-5): k := 0:
eq1 := diff(X(t), t) = -(A+B*X(t))*X(t):
ic[1] := X(365*k) = 1000:
s[1] := dsolve({eq1, ic[1]}, X(t), range = 0 .. .365, numeric);
pt[1] := eval([t, X(t)], s[1](0));
for k to 3 do
     tk := 365*k;
     A := rhs(s[k](tk)[2]);
     ic[k+1] := X(tk) = 500.*A;
     s[k+1] := dsolve({eq1, ic[k+1]}, X(t), range = tk .. 2*tk, numeric);
     pt[k+1] := eval([t, X(t)], s[k+1](tk));
end do;

plot([seq(pt[k], k= 1..4)], axes= boxed);

Use a table to store the solutions. Here's an example. I have a system of five randomly generated equations in variables v, w, x, y, z with a parameter a. I use a table indexed by a to store the fsolve solutions. Then I plot how x varies with a. (This particular plot is not very meaningful because I did nothing to control which of the system's many solutions was chosen by fsolve.)

eq||(1..5):= 'randpoly([v,w,x,y,z], degree=2)' $ 5:
for a from 0 by .01 to 1 do
     Sols[a]:= fsolve({eq1+a, eq||(2..5)})
end do:
plot([seq([a, eval(x, Sols[a])], a= 0..1, .01)]);

If you need more details, then please post your code.

The last command in the file is

slope[i][j]:= slope(pa[i], pa[j]);

In this command, you are using slope in two different ways---both as a procedure and as a place to store that procedure's results. This causes an infinite loop that eventually causes the kernel to crash.  You need to use a different name on the left side. You could make it SLOPE. Additionally, you may ultimately find it more convenient to use SLOPE[i,j] rather than SLOPE[i][j]. But I can't tell for sure without seeing what you intend to do with SLOPE.

Use RandomTools:-Generate and an `if` expression.

Tosser:= RandomTools:-Generate(float(range= 0..1, method= uniform, digits= 5), makeproc):
Toss:= ()-> `if`(Tosser() < .495, 0, 1):

Two examples of its use:
seq(Toss(), k= 1..10);
                  0, 1, 1, 0, 1, 1, 0, 1, 0, 1
add(Toss(), k= 1..10^5);
                             50325

Yes, the new ?LinearAlgebra,Generic package works over arbitrary commutative rings (many commands in the package require a field). There is no explicit rank command, but it's easy enough to get it from the row echelon form, for which there is a command. There's also SmithForm and BareissAlgorithm.

The help page ?LinearAlgebra,Generic shows an example of matrix computations over GF(3^3).

When dsolve (or any other Maple command) needs to make up a variable name which is not supplied by the user, it begins that name with _. These variables can be used like any other. The only significance of the _ is that it shows that the variable name was added by Maple.

You asked:

Does it mean that the integral can't be evaluated or maybe something else?

There are two separate issues here: the appearance of the _ in the solution and the appearance of an integral in the solution. As to whether the integral can be evaluated: It depends on the color of the integral sign. If it is blue then it can't be evaluated. If it is black, then it may be possible to evaluate it, and you can try to evaluate it using the value command. But, you don't need to remember the color code because you can apply the value command in either case.

DEtools[reduce_order](
     x*(x-1)*diff(y(x),x$2) - x*diff(y(x),x) + y(x) = 0,
     y(x)=x,
     output= solution
);

collect(expand(value(%)), [_C||(1..2)]);

You need to select based on type name rather than type symbol because indexed names such as a[1] are considered names but not symbols.

V:= < seq(`<|>`(selectremove(hastype, eq||k, name)), k= 1..3) >:
V1:= V[..,[1]]:  V2:= V[..,[2]]:
eq:= < < eq1 >, < eq2 >, < eq3 > >:

You need to use LinearAlgebra:-Equal instead of is, because is will not work with Matrices and Vectors.

LinearAlgebra:-Equal(V1+V2, eq);
                              true

You need to take a limit instead of using direct evaluation. Change the last line from

omega[L]:= 9:  t:= 10:  evalf(U[init]);

to the single command

limit(U[init], omega[L]= omega[G]);

and you will get 0.

Another option is to ?overload igcd so that it will take container objects (sets, lists, vectors, etc.) as arguments.

restart:
~set:= x-> convert(x, set):
local igcd:= overload([
     proc(X::~set, $) option overload; :-igcd(X[]) end proc,
     :-igcd
]):

igcd~([ [12,4], {8,2}, <35|5> ]);
                           [4, 2, 5]

See ?overload , ?coercion , ?updates,Maple17,LocalNames , ?dollar,Procedure (section "End-of-parameters marker $").