Do

remove(s-> true in map(e-> (is(rhs(e) < 0) assuming positive), s), Sa1);

Note for the experts: ormap won't work above because it refuses to handle the FAILs returned by is.

Short explanation: Only the right side of each equation is examined because it's assumed that the left sides are all simple names, as would be the case if these sets were returned by solve. For each right side, the is command tries to prove (using the ordinary rules of deductive logic and elementary algebra that we all know) that the expression must be negative under the assumption that all its names represent positive reals. If that proof succeeds (as indicated by returning true) for any case, the corresponding set containing that equation is removed. As is well-known and oft-discussed, is's ability to prove such things is quite limited, and it often returns FAIL, indicating that it can neither prove nor disprove a proposition. If this paragraph indicates that I've assumed something that you didn't intend, please let me know.

Use mul instead of product and add instead of sum because these operators allow the indices k to come from sets instead of ranges. The index sets are easy to construct with $ and minus. Like this:

restart:
N:= 4:
A:= Matrix(
   N$2,
   (i,j)-> if i=j then
      add(1/(x[i] - x[k]), k= {$1..i-1, $i+1..N})
   else
      mul(x[i]-x[k], k= {$1..N} minus {i,j})/mul(x[j]-x[k], k= {$1..j-1, $j+1..N})
   end if
);

If you're using 2D Input, the if ... end if might need to be changed to `if`(...) like this:

A:= Matrix(
   N$2,
   (i,j)-> `if`(
      i=j,
      add(1/(x[i] - x[k]), k= {$1..i-1, $i+1..N}),
      mul(x[i]-x[k], k= {$1..N} minus {i,j})/mul(x[j]-x[k], k= {$1..j-1, $j+1..N})
   )
);

(I'm not sure because I never use 2D Input.) Regardless, this latter form will work in any of Maple's interfaces.

Here is a procedure for it.

GramSchmidt:= proc(B::{Vector, list}, IP)
local n:= numelems(B), R:= Vector(n), M:= Vector(n), i, j;
   for i to n do
      R[i]:= B[i] - add(IP(R[j],B[i],args[3..])/M[j]*R[j], j= 1..i-1);
      if R[i] = 0 then error "Linear dependence detected at vector", i end if;
      M[i]:= IP(R[i]$2,args[3..])
   end do;
   <R|M>
end proc:
      
GramSchmidt(
   [seq(x^k, k= 0..5)],
   ((f,g,x::name)-> int((1-x)^(3/2)*f*g, x= 0..1)),
   x
);

The output is in two columns. The first is your orthogonal polynomials. The second is the constants that you called h_n.

You have dsolve({f1, bcs}, numeric), but it should be dsolve({f1, bcs[]}, numeric) because bcs is already a set. Also, you have a fourth derivative, so you need one more boundary condition.

The possibility of "getting an expression for each of them" seems to me exceedingly remote when you have those fractional exponents 0.62 and -0.05 on the dependent functions themselves. We need to go for a numeric solution, which can be obtained in seconds.

I need some clarifications:

• Is T supposed to be in the numerator or denominator of the exponent? The way that you wrote it puts it in the numerator. If it stays in the numerator, we need to rescale the units to maintain numeric accuracy because that makes it a negative exponent of very large magnitude. Putting T in the denominator, I got a smooth solution.
• You have the two derivatives equal to the same thing. Two functions whose derivatives are equal can differ by at most a constant. And since the initial conditions are equal, that constant will be 0. Is that what you intended? If so, there's no point in solving both as the solutions are guaranteed to be the same.

Some syntax corrections:

• x^-0.05 (for any expression x) needs parentheses: x^(-0.05).
• exp(1)^(x) (for any expression x) would be better written exp(x).

There is apparently a bug in HeatMap that prevents it from working when all the matrix values are the same. Dharr's workaround can deal with that. You can also do what you want by merging two sparsematrixplots, one which plots the 1s and  one which plots the -1s.

MyMatrixPlot:= (M::Matrix)->
   plots:-display(
      plots:-sparsematrixplot(M -~ (-1), color= blue),
      plots:-sparsematrixplot(M -~ 1, color= green)
   )
: 

Usage:

MyMatrixPlot(RM);

You asked

• Why is Maple definition of LegendreQ different from Mathematica definition?

This is explained is detail on the help page ?LegendreQ in the paragraphs beginning

"In standard references, the branch cuts are taken to be -infinity.. -1 and -1..1, and the values of the functions on the branch cut -1..1 are obtained by an averaging process (in order to obtain functions which are real valued on this interval).

"Rather than take this approach, Maple defines the default Legendre functions to be continuous from above onto the cuts -infinity..-1 and -1..1. However, alternative definitions of the Legendre functions are also available...." [emphasis added]

• And which is the correct one?

Let me ask you a philosophical question: Is it necessarily true that one is correct and the other is not?

• Also notice that the solution to this ODE should be real in |x|<1. So I do not undestand getting a complex solution.

How can you say that the solution that you got is not real without knowing the values of _C1 and _C2? If you specify any real initial conditions in the interval -1..1, these constants will be set so that the solution is real.

You also mentioned postfix notation. The output of InertForm:-Parse can be put into postfix notation with this procedure:

#Only intended to be used on InertForm output 
Postfix:= (f::function)->
   sprintf(
      "%A",
      evalindets(
         f, function, 
         f-> cat(nprintf("(%Q)", op(f)), nprintf("%A", substring(op(0,f), 2..)))
      )
   )
:

Your example:

Postfix(InertForm:-Parse("x+(y^(2+x)-4)/3"));

         "(x, (((y, (2, x)+)^, -4)+, 3)/)+"

Unlike the prefix output of InertForm:-Parse, the string output of Postfix has no inherent meaning to Maple, and it would take significant effort to write a procedure to automatically evaluate it.

As in the prefix form, the subtraction is handled as the addition of -4. The -4 itself is an integral operand; it is not the operand 4 acted on by an unary negation operator.

Here's one way with a single plot command. You can make many fine tunings of the ranges to see exactly the parts that you want.

eps:= .2:
x0:= Pi/4:
tan_line:= proc(f,x0) local x; unapply(D(f)(x0)*(x-x0)+f(x0), x) end proc:
plot(
   [tan, tan_line(tan, x0), arctan, tan_line(arctan, tan(x0))], 
   0..Pi/2-eps, 
   scaling= constrained, thickness= [2,0,2,0], view= [DEFAULT, -1..2]
);

We'll maintain two markers for the list, which I call lo and hi. Think of them as doing the same thing that your left hand does when you look up something in a paper dictionary: They keep track of the part that has yet to be searched. Then at each step, we examine the middle entry between lo and hi. If it's what we're looking for, then we're done. Otherwise, we can move one of the markers to the middle, thereby reducing the part that remains to be searched by a factor of a half.

The code below assumes 1) that the input list L is sorted with respect to the standard less-than operator `<`, 2) that the target x can be meaningfully compared via `<` to any entry of L, and 3) that if x is an entry, then it's an identical match as can determined by simple `=`. It'll be quite easy later to relax these restrictions.

BinarySearch:= proc(L::list, x)
local hi:= nops(L), lo:= 1, mid:= iquo(hi+lo, 2);
   while lo < hi do
      if x = L[mid] then return mid end if;
      if x < L[mid] then hi:= mid-1 else lo:= mid+1 end if;
      mid:= iquo(hi+lo, 2)
   end do;
   if mid <> 0 and x = L[mid] then mid end if
end proc:

Note that the second-to-last line could be written

if mid <> 0 and x = L[mid] then mid else NULL end if

However, there's never a need to write else NULL in Maple (or I believe in other Algol-based languages). If an if statement has no else clause, then it's assumed that it's else NULL.

Something like the above procedure is already implemented as ListTools:-BinarySearch. I thought that it'd be instructive for you to see the "bare bones" of this fundamental and very important algorithm.

Change every variable of the form A[b] to A__b. Note that there are two underscores. They will still display as subscripted. Using A[b] is unreliable because the meaning of this expression is not independent of A and b. On the other hand, A__b is totally independent from A and b.

If the solve still doesn't work after you do this, repost the newly subscripted expressions as a Reply to this Answer.

There are many ways to do those things. The commands that I show here are not the shortest ways, but I think that they are the best ways to learn first:

S:= seq(k, k= -10..20);
[S];
S:= seq(k, k= 115..231, 2);
[S];

The command unapply is the most general way to turn an expression into a function:

f:= unapply(Statistics:-PolynomialFit(4, X, Y, x), x);

Defining a function via

f(x):= ...

is unreliable, and clunky at best.

If an expression E contains the variable x, and E evaluates to itself, then you can do

f:= x-> E;

You must use exact arithmetic, not decimals, like this:

restart:
with(geom3d):
line(l1, [point(p1, [15, 6, 34/10]), vector([-4, 12, 3/10])]):
line(l2, [point(p2, [-17, 54, 32/10]), vector([6, -6, 2/10])]):
intersection(p, l1, l2):
coordinates(p);

The following formula is usually the first general formula for infinite summation that most of us learn, so it's likely that you're familiar with it:

sum(a*r^n, n= 0..infinity) assuming abs(r) < 1;

      -a/(r - 1)

Your first sum is clearly equivalent to sum(4^n, n= 1..infinity). If we ignore the abs(r) < 1 part and apply the above formula anyway, we get -4/(4 - 1) = -4/3. This is called "formal" summation. I don't understand why anyone would want this result, but I guess that some mathematicians have their reasons for wanting it since it has a name. Here's what the help page ?evalf,sum has to say about it:

"In the case of infinite sums, Levin's u-transform is used, which has the additional effect that sums that formally diverge may return a result that can be interpreted as evaluation of the analytic extension of the series for the sum (see the examples below).

"This behavior can be controlled through the use of the formal option or the _EnvFormal environment variable. If this variable is assigned the value false, or formal=false is specified, then this command will apply some (numeric) convergence tests to determine if the infinite sum in question is convergent, divergent, or otherwise non-convergent."

For your second sum, Maple doesn't know any way to simplify floor(5^n/4).