@nm Edgardo used the flat option to subsindets, and I did not. That's the cause of the difference you're seeing. Flat just means not recursive, which is the default. You'll need to decide which you prefer. I don't think that using the flat option is completely correct in this case: It's only performing your desired simplification to the outermost structures of the appropriate type. Edgardo should clarify his reason for doing this.

@Carl Love In the case you present in the PDF, all the items that are transformed are products of exp and sqrt (and possibly other factors). Thus they are all of type `*`. Thus we don't need to consider "naked" sqrts, those that aren't part of a `*` structure. Thus we don't need the _O or _1 factors. The procedure code immediately above simplifies to

MySimp:= expr-> 
    local rad:= anything^fraction;
    subsindets(
        expr, 
        And(
            `*`, 
            satisfies(
                e-> andmap(
                    membertype, 
                    [rad, {specop(exp), specop(exp)^anything}], 
                    {op}(e)
                )
            )       
        ),
        e->
            local d:= ilcm(map2(op, [2,2], indets(e, rad))[]);
            simplify((e^d)^(1/d))
    )
:

@Carl Love Your idea of bringing exp coefficients "under the radical" can be easily extended to radicals of arbitrary index:

MySimp:= expr-> 
    local _1, rad:= anything^fraction;
    eval(
        subsindets(
            subsindets(expr, rad, e-> _1*e), 
            And(
                `*`, 
                satisfies(
                    e-> andmap(
                        membertype, 
                        [rad, {specop(exp), specop(exp)^anything}], 
                        {op}(e)
                    )
                )       
            ),
            e->
                local d:= ilcm(map2(op, [2,2], indets(e, rad))[]);
                simplify((e^d)^(1/d))
        ),
        _1= 1
    )
:

@Carl Love Please let me know if you understand what I wrote above. If there are certain parts that you don't understand, let me know what those parts are. My explanation assumes some familiarity with "evaluation at a point" and its standard textbook notation; and some familiarity with the chain rule. I can back up the explanations to any particular point in your calculus studies. 

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