By trial and error, I found values of Digits and abserr to avoid those error messages. FWIW, those values were Digits=110 and abserr=100. This was only an initial exploration; I'm not suggesting that these would be good values to use in the finished product.

But, then a nonreal value was encountered. Looking at the ODEs, it's easy to see how that could arise from square roots and ln applied to the unknown function. There is no theoretical problem with a complex-valued BVP, but Maple cannot handle them directly. It can handle complex IVPs, which suggests that shooting may be a viable alternative. If it's tried, I suspect that the precision issues will become even more extreme.

If you believe (from theoretical considerations specific to your problem that I'm totally unaware of) that the expressions raised to fractional powers or that are arguments of ln in your ODEs should always be positive, then why that is not true needs to be investigated. It may be a result of round-off error and lack of precision.

@acer It's good to know that optimize treats a Matrix differently from a listlist. I wasn't suggesting that optimize not be used, only that the warning be ignored for the particular example shown by the OP.

Your instructor has asked you some fairly sophisticated questions about Maple's evaluation rules. Given that sophistication, I am dismayed and shocked that they think that it's acceptable to use the command matrix or the deprecated package linalg. They are mistaken; it's not acceptable. It's also disappointing that they give the impression that it's necessary to "load" a package before using it.

A lot of information that'll help you solve this can be found on help page ?eval.

A snag that you'll encounter if you use matrix​​​​​​ (instead of Matrix) is that a matrix uses the evaluation rule described on page ?last_name_eval. There's very little practical value in your instructor forcing you to learn that best-forgotten minutia of a long-dead command.

Unfortunately, whatever you learn in doing this exercise will be different if your variables are locals in a procedure. However, there is practical value in learning that.

@mmcdara In both of your failed color examples, you need to use ToPlotColor to convert the Color object into something understandable to plot commands. In your second example, with the named colors, you also need to use ColorTools:-Color.

I think that Maple 2019 is required for my code.

In the statement

for k,x in X do ...

k is the index of x in X, i.e., x = X[k] for every loop iteration. I think this was introduced in Maple 2019.

@adel-00 If you can't get a outside the integral, there's not much that can be done.

You wrote:

• The terms W1 and W2 are constant and assume precise values as a function of h.

Although I understood what you meant, at face value the above statement is self contradictory. A better wording would be "W1 and W2 are piecewise-constant functions of h." This is standard mathematical language; it is not derived from Maple's command piecewise.

@Spirithaunter You give up much too easily. There are a vast number of ways to do what you want. In particular, there is an error command, which will either stop the program or can be trapped with the try command to do anything you want.

MyCheck:= proc(F::{list, tabular}(realcons), UL::realcons)
    try
        if ormap(i-> is(F[i-1]+F[i] >= UL), [$2..numelems(F)]) then
            error "pairwise sum limit exceeded"
        else
            true
        fi
    catch:
        error 
            "expected 1st argument to be 1-dimensional with "
            "indices starting at 1"
    end try
end proc:

#Called via
MyCheck(F, 20);

You could also force a bad F to produce an invalid input error, like this:

TypeTools:-AddType(
    PairwiseSum,
    proc(F, UL::realcons)
        try
            F::{list, tabular}(realcons) and 
            andmap(i-> is(F[i-1]+F[i] < UL), [$2..numelems(F)])
        catch:
            false
        end try
    end proc
):
MyProc:= proc(F::PairwiseSum(20))
    #...whatever
end proc
:
F[1]:= 3: F[2]:= 20:
MyProc(F);
Error, invalid input: MyProc expects its 1st argument, F,
 to be of type PairwiseSum(20), but received F

@mmcdara I think that still only numeric edge weights are allowed. I recall writing an Answer here in the past few months with a kludgy workaround for that: You make the edge weights unusual numbers, numbers that won't otherwise appear in the problem. After calling DrawGraph​​​​​, substitute (using subs).the strings that you actually want for the string forms of the numbers. Of course, it works better when the old and new strings have the same length.

Although it's not realistic for an actual house, this problem supposes that the air doesn't move between rooms (convection), only the heat moves (conduction through the walls). So Newton's law applies to each pair of adjacent rooms (including the exterior), and the effect is additive.

@mmcdara The only thing that stopped me from answering this Question was that I didn't know how to account for the furnace F. Your solution says that the furnace raises the temperature of its surroundings at a constant rate, regardless of what that temperature is. My first guess was that it pumped out heat energy at a constant rate. Those two concepts are close to equal as long as the variation in the temperature of the surroundings isn't too extreme. And your way makes the ODEs linear, so it's likely what the professor intended. If you read the "boilerplate" on an actual furnace, the units will be energy/time (BTUs/hour in the US, I guess watts for the rest of the world) with an efficiency percentage given for the heat exchanger.

Anyway, vote up.

@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.