There's nothing wrong with your Maple: It can do the integral fairly quickly if you tell it, essentially, that the expression under the square root is nonnegative over the interval of integration. The conditions d/L__1 >= 2, d > 0 are sufficient. Use an assuming clause, like this:

simplify(...) assuming d/L__1 >= 2, d > 0;

where ... is exactly what you already have in the parentheses.

If these assumptions are not actually true for your problem, some deeper thought will be required.

It's nearly impossible to accurately compute the eigenvalues of a large matrix by finding the roots of the characteristic polynomial. Instead, use command LinearAlgebra:-Eigenvalues.

I show how to do that in this Post from a year ago: "Numerically solving BVPs that have many parameters" The BVP in that Post is extremely similar to yours; indeed, yours may be a direct simplififation of it. However, I wouldn't recommend that you attempt this if this is your first experience with Maple.

Like almost all operators in Maple, square brackets [] can be overloaded to have a special meaning for specific types of operands. Here is a simple implementation of what you're asking for:

Boolean_01:= module()
option package;
export `[]`:= proc(e::boolean, $) option overload; `if`(e,1,0) end proc;
end module:

with(Boolean_01);

#examples:
[3=7];
[3<7];

This is a minimalist implementation. There's no limit to the amount of subtlety that can be included.

1 /~ [8, 9, 9, 7, 9, 10, 5] mod 11

In general, to perforn an elementwise operation, append ~ to the operator. See help page ?elementwise.

The option scaling= constrained is the antithesis of being able to change the size of the z-axis. You should remove it.

Each seq command can only have one index variable. So, for example, if you wanted to generate all pairs of positive integers less than 4, you could use a double seq, like this:

seq(seq([i,j], i= 1..3), j= 1..3)

You could just as well iterate the j in the inside seq and i in the outside, which would generate the same pairs in a different order.

But you CAN'T do this:

seq([i,j], i= 1..3, j= 1..3)

simply because the syntax won't allow it, even though the concept is sensible. You attempted to do something like this, which was the reason for your error.

This is what I would do:

U:= eval(Array(1..2, 1..2, 1..2, 1..4, ()-> u[args]), Sol);

Now, U[2, 1, 1, 3], for example, would be the Array entry whose value corresponds to u[2, 1, 1, 3] in Sol. Unlike the other solutions based on either the assignment operator or assign command, this method does not cause any change to Sol or to the u's. Thus the u's can be used again as variables, perhaps to obtain another solution. This re-use of the u's would not affect the values now stored in the Array.

My method will work even if there are some index combinations for which there is no corresponding u in Sol.

My general policy is to avoid making any direct assignments to the fundamental variables of a problem (such as the independent and dependent variables); rather, I store solution values somewhere else, such as the Array, or just keep them in Sol for use on an as-needed basis. The command eval(expr, Sol) can be used to apply the values to any expression expr containing any of the u's.

In my code that you refer to, there is a line

Jxy:= ....

The purpose of that line is to list (in jet notation) the derivatives that you want to consider. In that case, I wrote some fancy code to extract all second derivatives of u and v with respect to x and/or y, but I could've just as well have written simply

Jxy:= {u[x,x], u[x,y], u[y,y], v[x,x], v[x,y], v[y,y]}.

Likewise, you can put whatever you want in the set. This is the only change that's needed.

Note that if the equation is nonlinear with respect to the desired derivatives, then a concept such as "the coefficient of u[x,x]" is not very meaningful. Nonlinearity of the overall equation is irrelevant; it's the linearity with respect to the desired derivatives that matters.

The problem is not recursive functions but a recursive assignment statement. Your first assignment statement inside the loop has x[n+1] on both sides of the assignment operator.

Assuming that I'm correct about the corrections that need to be made to the algorithm, as noted in my Reply above, then a translation of the algorithm into Maple is

Paar:= proc(M::Matrix)
local H, R:= copy(M), maxij, newcol;
     do
          maxij:= [max[index]((H:= rtable(antisymmetric, R^+.R)))];
          if H[maxij[]] <= 1 then return R fi;
          R:= <R|(newcol:= R[..,maxij[1]] *~ R[..,maxij[2]])>;
          R[..,maxij]:= R[..,maxij] *~  (1 -~ <newcol|newcol>)
     od
end proc:

 Note that I chose to return a new matrix rather than modify the input matrix M.

The above uses Maple 2018 syntax for embedded assignments. If you need to use it in an older Maple, let me know. It'll be easy to retrofit.

Use

op~([[1,2],[-2,1]], simplify(piecewise(op(f(t)), undefined)));
                                           [-1, 1]

We essentially have a pair of recurrences, one for a(n) and one for S(n). If we can eliminate one of two variables, then we can use the makeproc option of rsolve to get a numerical solver for the other. This is easy because a(n) = S(n) - S(n-1). Thus

SP:= rsolve({S(n)-S(n-1) = 2*S(n)^2/(2*S(n)-1), S(1)=1}, {S(n)}, makeproc):
SP(100)-SP(99);
                                          -2/39203

Note that no summation command is needed; there's not even one hidden in the procedure SP. I hope that Acer considers this to be the clever solution that he was expecting. 

The recurrence is nonlinear, so I don't have much hope for a closed-form solution.

You can extract the coefficients like this:

Jets:= {u,v}(x,y);
wave:= #You need to apply ToJet to the whole equation.
   ToJet(
      x*y*(diff(u(x,y),y,y) - diff(c(x,y)^2*diff(u(x,y),x),x)) - 
         TotalDiff(F(x, y, u[], u[1], u[2], v[], v[1], v[2]), y) + 
         TotalDiff(P(x, y, u[], u[1], u[2], v[], v[1], v[2]), x) = 
         0, 
      Jets
   );
#These are the 2nd derivatives that you want:
Jxy:= indets(wave, index~(op~(0,Jets), identical(op~(Jets)[])$2));
#The main coefficient-extraction command:
Cfs:= coeffs(collect((lhs-rhs)(wave), Jxy, distributed), Jxy, 'terms'):
#It needs to have its output put into a user-friendly form:
Cfs:= table([terms]=~[Cfs]):
#Now all coefficients are available by indexing Cfs:
Cfs[1]; #"constants" term
Cfs[u[x,x]]; #coeff of u[x,x] term

Where to start: Look up command ?addcoords