@Aschemer I am very curious about what is the optimal chunk size of the chunks from Y to use for this. Are the elements of Y generated iteratively? Are they stored on disk? It's clear to me that the automatic set sorting of X is beneficial for this. But it's not clear about the sorting of the chunks.

@raj2018 Any Reply posted to a Question's thread will put the thread at the head of the Recents queue, and thus anyone who pays attention to that will notice that there's a new Reply.

@student_md Looking at the inset integral formulas in your Question, they all essentially say this: "The symbolic integral of a product of 3 matrices and vectors equals that exact same product with a matrix of numeric constants thrown in the middle. The matrix of constants could be either the 2nd or the 3rd factor in the product." That doesn't make sense to me.

What are P(1) and P(2)?

@jrive I think that you meant to say that you disliked the position dependency of op, not that of indets, which is precisely the "master command" (within a family of commands) that allows you to avoid position dependency.

Anyway,, again I'm impressed that at this early stage in your Maple career you already realize some of the potential problems of position dependency. You may not realize this yet: The positions of many things can change from session to session. This makes position dependency much more of a problem than you otherwise realized. 

My most-common reaction to a parametric solution is "I don't want to see that much detail, although I do appreciate that it's there, and I certainly don't want to see a piecewise expression that's many screens wide." Filtering the output is tough work though, and it'd be great to have a subpackage to SolveTools for that. 

Now, I hope no-one takes that comment the wrong way, because the mathematical computation behind this is absolutely astounding, and the Maple code appears well written, although I've just skimmed it. It's in package RegularChains, which contains thousands of procedures (when you include the locals).

@nm I'm so sorry, change name to algebraic. I've been coding directly from my mind to my phone, with no Maple in between, the past few days because I'm on a long-distance household-moving project and getting very little sleep to boot. 

If you're syntax checking initial conditions from user input, you should check that they don't contain the independent variable (which would be nonsense IMO). Thinking about that is why I had name on my mind.

If you display your worksheet inline here, there's a good chance that I could answer your question immediately from my phone. I'm traveling at the moment, and so I can't load your linked worksheet into Maple to view it.

@emendes I have no question about either the syntax or semantics of @, so there's not much point in you making any major accuracy test of it. For any f, g, and x (even if they're sequences), (f@g)(x) = f(g(x)) and (f@g)~(x) = (f~@g~)(x). (In both cases, that's identical equality, not just mathematical equality.) The only question that I have is whether either of the forms on the left sides of those two are more efficient than the corresponding right sides, and whether that may depend on whether f or g are built-in. 

Regarding the rational functions: Yes, that's exactly what I meant.

@Carl Love My Reply immediately above has been significantly updated, with very important details. I only added stuff without altering its original content. So, I'm just posting this to bump this thread up in the Recent queue.

@vv I'm speaking of what may be possible mathematically, not of what can be done with a stock Maple command. If one knows only a finite number of terms of two series f and g, there are many functional operators F such that one can find a finite number of terms of F(f,g). I think that something akin to that may be possible. 

@emendes Use this as a replacement for abs:

Pabs:= e-> {e, -e}[1]:

Then, regardless of whether you use Classify or my procedure, replace rhs~ with (Pabs@rhs)~.

Given two single-argument commands f and g, I'd guess that (f@g)~(e) is more efficient than f~(g~(e)), and I'm even more inclined to make that guess if at least one of f or g is built-in (as is rhs). The reason that I'd guess that is that I think that `@` itself is built-in, giving it an ability to optimize this in ways that are not possible for Maple-level code. I'm traveling and on my phone right now, so I can't check this, but @acer or @Joe Rlel may know off the top of their heads.

Both these or the versions with simply rhs will only detect equalities that can be detected immediately, i.e., without doing any symbolic work such as simplify. Since I'd generally expect that different solutions to a particular system of equations would be returned in the same general form by solve, I don't expect this to be a problem. (For example, it seems unlikely that 2^(1/2) and 4^(1/4) would both appear in the same solution set.) Are the initial equations always rational functions with rational coefficients? If so, then it seems even more unlikely.

@Maple_lover1 Yes, here is my version that makes those requested changes, as well as a few others. I use the same subs-into-a-procedure-template technique that acer used in his Reply "Another way". This technique is extremely powerful symbolically (much more than unapply), yet also extremely efficient; however, it takes a lot of practice to master.

By replacing the doubly nested procedures with triply nested procedures, I was able to isolate the second simplify so that it would be performed only once for each n and m, but with symbolic x. I also was able to replace the "heavy" command piecewise with a simple `if`.

functionB:= proc(
    k::nonnegint, M::nonnegint,
    epsilon::(positive &under (evalf@`+`, 1)), gamma::(positive &under (evalf@`+`, 1))
)
option remember;
local m;
    subs(
        [
            __K= 2^(k-1), __M= M, 
            __h= subs(
                m= __m,
                simplify(
                    2^(epsilon+gamma+1)*GAMMA(epsilon+m+1)*GAMMA(gamma+m+1)/
                        (2*m+1+epsilon+gamma)/m!/GAMMA(epsilon+gamma+1),
                    assume= m::nonnegint
                )
            )
        ],
        proc(n::integer[0..__K], __m::integer[0..__M-1])
        option remember;
        local x;
            subs(
                 __f= subs(
                     x= __x, 
                     simplify(2^(k/2)*__h*JacobiP(__m, epsilon, gamma, 2^k*x-2*n+1))
                 ),
                 (__x::algebraic)-> `if`(n-1 <= __x*__K and __x*__K <= n, __f, 0)
            )
        end proc
    ),  
    (x::algebraic)-> (1-x)^epsilon*(1+x)^gamma
end proc
:

I tested this identically to how acer tested his, including the initial assignments to globals n, m, and x. I obtained identical results, except for the wording of the error messages (which I strongly prefer to handle through argument type checking). Because of the triple nesting, all calls such as psi(-1, 0, x) need to be changed to psi(-1,0)(x).

All variables that begin with double-underscore are what acer called "dummies". They are the variables within the procedure templates which receive the substitutions from outside the templates. Unfortunately, the subs technique requires that these be global unassigned symbols. Since you should never make an assignment to any global variable beginning with a double underscore, this is likely safe. When the subs technique is used correctly, none of these globals remain in the returned procedures.

To facilitate argument-based type checking, I used two parameterized structured types that are unfamilar to most readers. The first is &under, and the second is integer[A..B]. 

type(e, T &under f) is true iff type(f(e), T) is true.
type(e, T &under (f, a)) is true iff type(f(e,a), T) is true.

type(n, integer[A..B]) is true iff n::integer and A <= n and n <= B.

@Kitonum You wrote:

• Since Maple does not explicitly solve your system of differential equations (only numerically)....

That's not quite true. Maple gives an exact series solution, which is not as sophisticated as a fully symbolic solution, but more sophisticated than a numeric solution. There's a lot of symbolic processing that can be done with series. If we suppose the existence of a function F such that y(t) = F(x(t)) for all t (possibly restricted to some interval), it may be possible to determine the coefficients of series(F(x), x= ...) (i.e., with x as the independent variable). I don't know how to do this for this particular system of equations, but I'm not yet willing to say that it's impossible, unless you have some proof (or even evidence). 

@cbowers Try changing ^+ to ^%T. Both are operators for the matrix/vector transpose. The %T may work in your Maple version. If not, let me know, because there are other easy options.