@vv As is almost always true in Maple, the = 0 is not needed; (lhs-rhs)(mysol) is sufficient.

@T K Burch If you have x and y data in separate column vectors X and Y (and, as before, F is an expression with a single free variable, x), then you can do this:

plot([F, <X|Y>], x= (min..max)(X), style= [line, point])

Column is the default vector orientation. If for some reason you have row vectors, let me know. The <...|...> operator is a columnwise matrix constructor.

@T K Burch So, this small thing of the punctuation at the ends of statements makes a big difference in Flow. In Maple, the colon means to suppress the output. The semicolon means to not suppress it, and is not required unless there's an immediately following statement.

If Flow doesn't like the punctuation, it should give a normal error message about it, not the weird message that you reported. STATSEQ stands for "statement sequence", which indicates that it was expecting another statement after the punctuation, which is understandable. subsindets stands for "substitute the indeterminants". It's one of the most important low-level Maple commands, but it's a command for a programmer, not an end user.

So, it seems like the Maple Flow input is being processed by regular Maple code. That's useful to know. 

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