@nm As I said to you in another recent thread:

• This hasn't come up in this thread yet, but when doing this type of work, it's important to keep in mind that Maple's internal representation of a/b is a*b^(-1).

And, indeed, it didn't come up in that thread, but now it's an issue here.

In our common handwritten and typeset mathematical notations, fractions appear to be an operator with two operands, numerator and denominator. And, indeed, that is literally true when they're encoded in Latex or "markup" languages. But Maple doesn't store them that way (for good reasons).

So, you can correct this issue several ways, including those given by Edgardo and acer. My first suggestion is to replace specop(exp) with {specop(exp), specop(exp)^anything}.

Given types t1 and t2, the type {t1, t2} is almost the same as Or(t1, t2). The only difference is that with Or, you can guarantee the order that the types are checked. That's critically important when constructing recursive types (not something that we're doing here).

@John Eric Your final expression Rt contains two free variables: theta and h[1][1]. I suspect that the h[1][1] was unintended.

@Gabriel samaila Okay, you've made some (small) changes to the system of equations, which is one reason it wasn't working.

The other reason is that the || catenation operator doesn't work (at least in this case) with 2D Input, where it's interpretted as a vector norm operator. So, change this line

sys:= {bcs||(1..2), ode||(3..4)}:

to

sys:= {bcs1, bcs2, ode1, ode2}:

You changed fixedparameter:= [pr= 1] to [N= 1]. This makes no sense in the context of the plot that we want right now! We need to vary N (as Nu) along the horizontal axis. But pr needs a value, and I just used pr=1 again.

Look at the plot in my Answer above. Notice that the entire variation of f''(0) is in the 10th decimal place. The new plot is very similar, except that (if I did these computations to sufficient accuracy) the entire variation is in the 21st decimal place! It requires very high accuracy in the dsolve command to do this. I used

Digits:= 25:
and I added the options
maxmesh= 2^15, abserr= 1e-15
to the dsolve command.

I don't know if that is sufficient accuracy! I think maybe Digits should be higher (try 30) and abserr should be lower (try 1e-23). I don't think that you'll need to change maxmesh. These make generating the plot take an extremely long time.

I added a Reply to my Answer to cover this new situation.

Okay, based on your comments elsewhere, I think that you want to select the operands that contain a unless they also contain exactly one other variable. In that case, use

select(t-> has(t,a) and not nops(indets(t, name))=2, [op](_x+f));

@nm You're right that my version will remove a/b, a^b, etc. I assumed that that was what was wanted. But, the OP's specifications are very unclear.

This hasn't come up in this thread yet, but when doing this type of work, it's important to keep in mind that Maple's internal representation of a/b is a*b^(-1).

I'm excited to see these new variations. I hope that you post them here soon.

@nm Sorry, my mistake. Change the subsindets command to

subsindets(
    subsindets(expr, sq, e-> _O*e), 
    And(`*`, satisfies(e-> membertype(sq, {op}(e)))),
    e-> LargeExpressions:-Veil[Z](eval(e, _O= 1))
);

The trick is to make a preliminary alteration to the expression so that no sqrt is "naked", i.e., they're all of type anything*sqrt. Then in the outer subsindets we only look for that form, which'll now catch all the sqrts, even those that were originally naked.

I'll make this correction in the Answer also.

@nm The Mathematica command that you show is obviously equivalent to subsindets, which is a pattern-matching command; the Mathematica command is not equivalent to patmatch. I don't have anything against pattern matching in general; I just think that its Maple implementation in patmatch, applyrule, etc., is junk. Those commands are to my mind obsolete, and their design was half-baked in the first place.

I think you and others are attached to patmatch and applyrule simply because their names are more meaningful than subsindets, not because they have superlative functionality.

In my opinion, subsindets (or evalindets) or its equivalent in Mathematica or other languages is the most important command in all of computer algebra.

I thought that it'd be more generally useful to separate the cases as Z[1], ..., Z[5]. I don't consider ee->Z to be "the normal transformer" because I don't think that it's normal to replace distinct subexpressions with the same variable. This way, you can examine the contents of all the Zs via

LargeExpressions:-Unveil[Z]([seq](Z[k], k= 1..LargeExpressions:-LastUsed[Z]));

But separate Zs isn't really necessary. If you want to replace all the expressions with just Z, you could do that. However, note that that could lead to unexpected oversimplifications such as the cancelling of two Zs which originally represented distinct expressions.

@lcz It's not clear what you want. Why should a*b*c be kept? Do you only want to remove terms containing a and exactly one other variable?

@MapleEnthusiast Why do you care about an iterative solution? If you want a solution equivalent to the Moore-Penrose, then use LinearAlgebra:-LeastSquares (which'll do it without explicitly converting the pseudoinverse). If you want to know the spectral radius, just compute that.

AFAIK, an iterative solution to a linear system is only worthwhile if the coefficient matrix is large and sparse. The iterative process that has been described so far in this thread (which I guess could be called fixed-point iteration) is crude and not often used.

@Carl Love The labels can also be generated on an as-needed basis from their position number by an arithmetic formula:

GetLabel:= subs(
    _AZ= ["", cat("A".."Z")],
    proc(N::integer[1..2^14])
        local n:= N, k;
        cat(seq(_AZ[1+iquo(n,26^(2-k),'n')], k= 0..2), 1)
    end proc
):
GetLabel(1234);
             "AUL1"

In Maple 2019 or later and 1D input, the procedure can be shortened to

GetLabel:= subs(
    _AZ= ["", cat("A".."Z")],
    (N::integer[1..2^14])->
        local n:= N, k:= 3;
        cat('_AZ[1+iquo(n,26^--k,'n')]'$3, 1)
):

@Carl Love Recent versions of Excel allow up to 2^14 = 16,384 columns by using up to 3 letters for the labels. The code above can be easily extended to accomodate all of them:

my_alph:= [seq](cat(("A".."Z")$k, 1), k= 1..3)[..2^14]:

While 16000 may seem like a lot to us, the execution time of the above is barely even measurable by Maple.

@emendes In my opinion, alias is superficial. Yes, it can be used to achieve a certain desirable output display, and it allows for abbreviated input, but it doesn't fundamentally affect the computation. While the bugs are still being worked out, the presence of alias just adds confusion. Please present the above Lie derivative example from restart with any alias commands commented out (not removed). They can be replaced after the bugs are worked out.