Go to the Maple Applications Center to get the Orthogonal Expansions package. Click on "Download attached file". Follow the trivial installation instructions in the ReadMe file.

Now you can easily find Fourier series in Maple. (Note that Pi is spelled with a capital P in Maple.)

restart:
f:= x-> min(abs(x), Pi/2):
FS:= OrthogonalExpansions:-FourierSeries(f(x), x= -Pi..Pi, infinity);

plot(f(x), x= -Pi..Pi);

plot([seq(eval(FS, infinity= n),  n= [1,2,3,5,6,30])], x= -Pi..Pi);

You will gain understanding of the order in which things are done if you set printlevel:= 15 before running your loop. You are reusing k for every i and j. So only the last (i,j) combination is saved for each k. If you don't understand that, I can explain more. Let me know.

What you want is this:

for i to 3 do
     for j to 3 do
          k:= 3*(i-1)+j;
          a[k]:= (i,j)
     end do
end do;

eval(a);

Note that table entries are not necessarily stored in any order.

Your nested for-do loops can also be done by nested seq loops and a table constructor:

a:= table([seq(seq(3*(i-1)+j = (i,j), j= 1..3), i= 1..3)]);

Try

subs(sqrt(c)= omega*sqrt(m), xsol);
simplify(%);

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