@Preben Alsholm There are 32 unknowns and 24 equations. The number of solutions might depend on which 8 variables solve decides to use in the basis, and that may change on different runs.

I tried to solve for arbitrary matrix T as the Asker's updated post asks. I used eliminate. I let it run for 2 hours and use 7 Gigs memory, and then I killed it.

Try raising your cutoff value of 2.5e-2 to 0.3 or 0.4. Then you can really see the points from which it is interpolating the surface. You can see a collection of points on each side of the surface such that the surface lies neatly between them. It might also be instructive to vary the symbolsize of your plotted points, with points closer to the surface plotted larger. This might give more of a clue as to how the algorithm works.

The trick is to make each frame include a static (i.e. not insequence) display of all previous frames. This can be done with any animation. Here's the original code (which I indented) and my modifications. I added a splash of color. I also changed the sum to add for efficiency, which necessitated making g a procedure.

(***** begin original (comment out) *****
g:= sum(arctan(1/sqrt(n)), n = 1 .. i):
p:= seq(
    plots[polygonplot](
        [[0, 0],
        [sqrt(k+1)*cos(subs(i = k, g)), sqrt(k+1)*sin(subs(i = k, g))],
        [sqrt(k)*cos(subs(i = k-1, g)), sqrt(k)*sin(subs(i = k-1, g))]
        ], color = white
    ), k = 1 .. 70
):
plots[display]([p], insequence = true);
***** end original (comment out) *****)

N:= 70:
g:= k-> add(arctan(1/sqrt(n)), n= 1..k):
p:= [seq](
    plots:-polygonplot(
        [[0, 0], sqrt(k+1)*~[cos,sin](g(k)), sqrt(k)*~[cos,sin](g(k-1))    
        ], color= COLOUR(HUE, k/N)
    ), k= 1..N
):
plots:-display([seq](plots:-display(p[1..K]), K= 1..N), insequence = true);

The trick to using evalhf here is that the integration needs to be done at the time h is created, not at the time that h is called. And there's extra complexity that n = 2 needs to be handled as a separate case to avoid division by zero. So, you need to define h this way:

h := subs(J= int(f(x)*sin(n*x), x= -Pi..Pi)/Pi, n-> `if`(n=2, 1, J));

Your f and g can remain the same. If you define h as above, then the call evalhf(g(2)) will work. This is correct usage, not effective usage. To use evalhf effectively, you need to pack as much computation as possible into each invocation of evalhf.

I can't answer the "why?" of your question, but there is a very easy workaround: Make the order 2 or more.

Unfortunately, freeze doesn't work for integers and rationals. But subs does work. So the following kludge does it:

ex:= 16^j/16;
exf:= subs(_x= 16, simplify(subs([1/16= 1/_x, 16= _x], ex)));

It sure looks ugly though.

You have a line

H := dsolve*(ivp, {e(t), w(t)}, numeric);

The asterisk in that line should be removed.

Also, it is not an error, but there is no good reason to use `union` in prefix form. You might as well write

ivp:= sys union ic;

which is closer to the usual mathematical notation.

There are many ways to plot it. Here are two:

plots:-odeplot(H, [[t,e(t)],[t,w(t)]], t= 0..100, color= [blue,red]);

plots:-odeplot(H, [e(t),w(t)], t= 0..100);

Why do you want it to be an array of functions? Why not just have phi be a function of two variables, i and x?

phi:= (i,x)-> (i+2)/(i+1)*x^i + x^(i+1);

If you want to get fancy with the syntax and still be able to invoke phi with the square brackets, you can do this:

phi:= proc(x)
   local i;
   i:= op(eval(procname));
   (i+2)/(i+1)*x^i+x^(i+1)
end proc:

phi[2](x);

Both of these approaches have the advantage of not polluting the memory with n copies of what is essentially the same function.

This works:

DE := [diff(x(t), t) = sin(x(t))/sin(t)];
with(DEtools):
DC:= sin(x(t))/sin(t):
dfieldplot(DE, x(t), t = -2 .. 2, x = -1 .. 2, arrows = medium, title = 'Diff', color= DC);

Note that x is always x(t) in the differential equation. What were you trying to achieve with the .1 in the color?

The command is solve(..., B). It is not clear what you mean by h and by "equation". Is h meant to be a variable in the equation? Then your := should be =. Or do you mean that the expression assigned to h should be set to 0? Either way, the answer is not so simple: You have a degree-6 polynomial in B.

Notation: I'll use juxtaposition to indicate intersection, `+` for union, `-` for minus, with `+` having precedence. |A| is nops(A), and (n C k) is a binomial coefficient.

We split cases via |AB|. There are three cases, |AB| = 3, 4, or 5. Regardless of the case, we have a factor of (30 C 5) for the choice of A+B and a factor of 2^25 choices for C - A+B.

Case |AB| = 5: Thus A=B. There are (5 C 3) = 10 choices for ABC in AB.

Case |AB| = 4: Either |A|=4 and |B|=5 or vice versa. This dichotomy contributes a factor of 2, and we will consider only the |A|=4 case without loss of generality. There are (5 C 4) choices of A within A+B and (4 C 3) choices of ABC in A. Total factor for this case 2*(5 C 4)*(4 C 3) = 40.

Case |AB| = 3:

• Subcase |A|=3, |B|=5 or vice versa. Factor of 2 for the dichotomy. (5 C 3) choices of A in A+B, and only 1 choice of ABC in A. Total factor: 20
• Subcase |A|=4, |B|=4. (5 C 4) choices of A in A+B, (4 C 3) choices of AB in A, and(3 C 3) choices of ABC in AB, for a total of 20.

Grand total: (20+20+40+10)*(30 C 5)*2^25 = 430353709793280 =

Edit: Corrections to the above.

Case |AB| = 4: Either |A|=4 and |B|=5 or vice versa. This dichotomy contributes a factor of 2, and we will consider only the |A|=4 case without loss of generality. There are (5 C 4) choices of A within A+B and (4 C 3) choices of ABC in A. The choice of whether the lone element of B-A is or is not in C contributes a factor of 2. Total factor for this case 2*(5 C 4)*(4 C 3)*2 = 80.

Case |AB| = 3:

• Subcase |A|=3, |B|=5 or vice versa. Factor of 2 for the dichotomy. (5 C 3) choices of A in A+B, and only 1 choice of ABC in A. For each of the 2 elements of B-A, there is a choice of whether it is in C, for a factor of 2^2. Total factor: 2*(5 C 3)*2^2 = 80.
• Subcase |A|=4, |B|=4. (5 C 4) choices of A in A+B, (4 C 3) choices of AB in A, and(3 C 3) choices of ABC in AB. For each of the 2 elements of (A-B)+(B-A), there is the choice of whether it is in C. Total factor: (5 C 4)*(4 C 3)*2^2 = 80.

Grand total: (10+80+80+80)*(30 C 5)*2^25 = =

There must have been a floating-point number introduced into the computation, either explicitly or through use of evalf or a similar command. What were the steps that led to your 2.*I*q*r? This is a phenomenon called "floating-point contagion", (which makes it sound like a bug, but it it the intended behavior). Basically, the decimal points spread from number to number in your expression. To get rid of them,

subsindets(ex, float, convert, rational);

where ex is the expression that you want to correct.

You can't use subs to substitute for something that doesn't exist as a distinct entity. Notice the difference in the following two uses of subs.

To be precise about what I mean by "distinct entity" requires some knowledge of Maple's internal representations of objects.

But, we can ignore that, because there is a command that is much more sophisticated than subs, and which doesn't care about the internal representations. That command is algsubs. It only does one substitution at a time. So, to use it for your case, we need to put it into a loop.

Replace the line

Eqs1:= subs(BC, Eqs);

with

for s in BC do Eqs:= algsubs(s, Eqs) end do; Eqs1:= Eqs;

I did this in your code, and it ran through to the end, generating the point plot.

Assuming that ex is set to the expression that you want to manipulate, it can be done as

sort(-numer(-ex)/denom(ex), x);

The second argument to sort, the x, is to tell it which variable that you want to appear first.

I can show you how to do it by TypeTools if that's what you really want. But isn't this preferable?

select(has, G1, w);

Let me know.

Also, you are using the wrong quotes again, you mean 'w', not `w`. The latter does nothing to something which is already a name. As for the '[]', I doubt that that's what you really want. There's `[]`, [], and [][]. They all do different things when applied as functions.