Please read carefully the section "Using Assumptions" in the help page ?solve,details.

Conclusion: Don't use assuming real with solve.

After what acer taught you yesterday, I'd've thought that you'd stop using all-variable assumptions such as assuming real for any command. 

This code uses the extremely fast package LinearAlgebra:-Modular. It is much faster than the compiled code, finishing the 2^25 case in under 2 minutes.

[Edit: The complete corrected code is in a Reply below. It does the 2^25 case in 0.2 to 0.3 seconds.]

BooleMobiusTransform:= proc(V::Vector(datatype= integer[8]))
uses LAM= LinearAlgebra:-Modular; 
local 
    nv:= numelems(V), n:= ilog2(nv), im:= Scale2(1, n-1), istep:= im, 
    jm:= 1, h:= Scale2(im, 1), istart, 
    R:= LAM:-Create(2, h, 0, integer[8]);
    LAM:-Copy(2, V, 1..-1, R, 1..nv);
    to n do 
        istart:= 1; 
        to jm do
            LAM:-AddMultiple(
                2, R, istart..im-1+istart, 
                R, istart+istep..im-1+istart+istep
            );   
            istart+= h
        od;
        im:= Scale2(im,-1); istep:= Scale2(istep,-1); 
        jm:= Scale2(jm,1); h:= Scale2(h,-1)
    od; 
    R
end proc
:
for n from 11 to 19 do
    V:= Vector(2^n, rand(0..1), datatype= integer[8]):
    lprint('n'=n);
    V1:= CodeTools:-Usage(BooleMobiusTransform(V))
od:
n = 11
memory used=0.58MiB, alloc change=0 bytes, 
cpu time=15.00ms, real time=9.00ms, gc time=0ns
n = 12
memory used=1.16MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=14.00ms, gc time=0ns
n = 13
memory used=2.35MiB, alloc change=0 bytes, 
cpu time=31.00ms, real time=21.00ms, gc time=0ns
n = 14
memory used=4.76MiB, alloc change=0 bytes, 
cpu time=47.00ms, real time=44.00ms, gc time=0ns
n = 15
memory used=9.63MiB, alloc change=0 bytes, 
cpu time=93.00ms, real time=85.00ms, gc time=0ns
n = 16
memory used=19.51MiB, alloc change=-1.00GiB,
cpu time=390.00ms, real time=358.00ms, gc time=93.75ms
n = 17
memory used=39.51MiB, alloc change=1.00MiB,
cpu time=437.00ms, real time=393.00ms, gc time=109.38ms
n = 18
memory used=80.01MiB, alloc change=2.00MiB, 
cpu time=797.00ms, real time=742.00ms, gc time=125.00ms
n = 19
memory used=162.01MiB, alloc change=3.00MiB, 
cpu time=1.62s, real time=1.50s, gc time=265.62ms

Often (but not always), Maple can symbolically integrate the Heaviside forms of piecewise expresssions that it can't do directly. Two branches of your proposed result needed small corrections (assuming I and Maple did everything correctly here).

restart:
interface(showassumed= 0):
assume(n::posint, j::nonnegint, j <= n, T > 0, alpha > 0);
h:= T/n:
psi:= t->
    local i:= op(procname);
    convert(
        if i=0 then piecewise(t<0, 0, t<h, 1-t/h)
        elif i=n then piecewise(t < T-h, 0, 1-(T-t)/h)
        else 
            piecewise(
                t < (i-1)*h, 0, 
                t < i*h,     t/h-(i-1), 
                t < (i+1)*h, (i+1)-t/h
            )
        fi,
        Heaviside
    )
:    
g:= proc()
local i,j,t;
    (i,j):= op(procname);
    simplify(
        int((j*h-t)^(alpha-1)*psi[i](t), t= 0..j*h)/GAMMA(alpha)
    )
end proc
:
g[0,j]() assuming additionally, j > 0;
             1         //       alpha            (alpha + 1)
    - ---------------- \\(j - 1)      (j - 1) + j           
      GAMMA(alpha + 2)                                      

                       alpha\  (-alpha)  alpha\
       + (-alpha - 1) j     / n         T     /


#Compare with your proposed result:
#Note that I added "-" in front of (j-1)!
simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        (-(j-1)^(alpha+1) + j^alpha*(1-j+alpha))
);
                               0

g[j,j]() assuming additionally, j > 0;
                         (-alpha)  alpha
                        n         T     
                        ----------------
                        GAMMA(alpha + 2)

#Clearly, that's your proposed result.

g[i,j]() assuming i::posint, additionally, i > j, 
    additionally, i <= n;
                               0

g[i,j]() assuming i::posint, additionally, i < j;
       1         / alpha /           alpha            
---------------- \T      \(j - i - 1)      (i - j + 1)
GAMMA(alpha + 2)                                      

                              alpha            alpha            \ 
   + (-i + j + 1) (-i + j + 1)      + 2 (j - i)      (alpha + 1)/ 

   (-alpha)\
  n        /


#Note how I corrected your formula below:

simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        ((j-i+1)^(alpha+1) + 2*(j-i)^alpha*(alpha+1) - (j-i-1)^(alpha+1))
);
                               0

It can be done like this (Warning: The key line of this code contains a large number of weird looking punctuation marks):

SpecialIndexPrint:= (e::uneval)->
    subsindets(
        convert(op(0, eval(e,1)), `local`)[eval([op](eval(e,1)))[]] = 
            eval(e),
        typeindex(symbol),
        n-> nprintf(
            cat(`#msub(mi(\"%a\"),mo(\"`, `%a ` $ nops(n)-1, `%a\"))`),
            op(0,n), op(n)
        )
    )
:
for i to 2 do  
    for j to 2 do 
        S[i, j] := 2*mu*varepsilon[i, j] - add((2*mu)/3*varepsilon[r, r]*delta[i, j], r = 1 .. 3); 
        print(SpecialIndexPrint(S[i, j])); 
    end do;
end do:

Yes, my procedure SpecialIndexPrint is more complicated than the Answers by Kitonum and Tom, but it does have a few advantages:

• It can be used either inside or outside the loop.
• It is not specific to the name S.
• It works on any indexed symbols (with any number of indices).
• It replaces the commas with spaces (although it'd be easy to add a correction for this to either of those other Answers).
• It replaces only the commas that occur as index separators.

The prettyprinted results of the above cannot be adequately displayed here in plaintext form. You'll need to execute it in a Standard worksheet to see it. But it's exactly what you asked for.

To better understand how the above procedure works, use lprint instead of print, for example,

lprint(SpecialIndexPrint(S[1,2]));

`#msub(mi("S"),mo("1 2"))` = 2*mu*`#msub(mi("varepsilon"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("1 1"))`*`#msub(mi("delta"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("2 2"))`*`#msub(mi("delta"),mo("1 2"))`-2/3*mu*
`#msub(mi("varepsilon"),mo("3 3"))`*`#msub(mi("delta"),mo("1 2"))`
 

Maple doesn't "simplify" equations, and from what you've shown of Mathematica (in this and in previous Questions), it doesn't either. When you use simplify on an equation, it just simplifies the two sides separately. Thus, the following can be and should be used for all equations, regardless of whether they have square roots:

(numer@simplify@(rhs-lhs))(new_ode);
                       3                          
             / d      \        / d      \         
             |--- u(x)|  + 8 x |--- u(x)| - 8 u(x)
             \ dx     /        \ dx     /         

In American algebra education, the word simplify is not applied to equations; one simplifies expressions (and equations are not considered expressions), and one solves equations. I assume that this is the same in Kitonum's country, so I agree with what he said about you confusing simplify and solve.

However, note that my solution above does not use solve! So I think that it is a black-box solution such as you need. Instead, the equation is converted to an expression, which we thus assume is equal to 0. The expression is simplified, and since it equals 0, we're justified in further simplifying it to just its numerator.

You'll make your coding much simpler if you immediately process all input equations with numer@(lhs-rhs). Indeed, in this case, the simplify isn't even needed: (numer@(lhs-rhs))(new_ode) produces the same output. I just included the simplify for more general cases. I would guess that this process of converting equations into expressions is nearly universal in computer algebra. Certainly it's older than Maple.

Using the type &inside that I defined in the previous Answer, all you need is this:

evalindets(expr, specfunc(abs) &inside specfunc(ln), eval, abs= _@@0)

_@@0 is a fewest-characters representation of the multivariate identity function; you could also use abs= (x-> x).

Targetting only 1-argument abs is possible, but a little more complicated.
 

You can create your own structured types, and a generalization of what you describe is something that I want often enough that it's worth having a "prepositional" structured type for it. Specifically, I want to check for instances of type A inside type B.

TypeTools:-AddType(`&inside`, (e, `in`::type, out::type)-> e::out and hastype(e, `in`)):        
expr:=sin(x)+ln(abs(x))+ln(x+1/sqrt(abs(x+3)))+ln(x^3):

indets(expr, specfunc(abs) &inside specfunc(ln));

It's been several years since anything has been required as a first argument of specfunc. The first argument defaults to anything if none is supplied.

I do appreciate the ambiguities that acer pointed out; however, having read and Answered several hundred of your thoroughly written^[*1] Questions (which I almost always Vote Up even if I don't Answer), I'm pretty sure that I understand what you mean:

• You know the specifc meaning of function as a Maple type (and due to a recent Answer you know that the evaluated result of sqrt is not a function).
• You want the recursive scan that indets (and its partner hastype) usually provide.

[*1] Although your English grammar, spelling, etc., could be improved, I still almost always understand what you mean. I do realize that English is probably not your native language. I'd be happy to give ad hoc grammar advice if that'd be appreciated.

Always use Matrix instead of matrix.

Please post your worksheet and state your Maple version.

I believe that you should use D, but only indirectly, called from your own procedure. Rather than loading D's remember table, my procedure uses a separate table that only holds your special derivatives. This handles both symbolic and numeric constants correctly. It's also a very short code. The fact that my table is declared sparse means that the derivatives of any symbols not specifically indexed (which are presumed to be constants) will be 0.

restart:
MyDerivatives:= table(sparse, [x= x*y, y= 1, z= z^2]):
MyD:= e-> evalindets(D(e), specfunc(D), d-> MyDerivatives[op(d)]):
  
MyD(x*y^2/2 + 3*y^3*z*b + 3 + a);
              1    3            2          3  2  
              - x y  + x y + 9 y  z b + 3 y  z  b
              2                                  

Here are the "under the hood" basics of how assuming works. Understanding this should clear up your confusion. It is obvious that the proper evaluation of F(G(e1, e2)) assuming a1, a2 requires that the assumptions be applied before F(G(e1, e2)) is evaluated. In other words, the normal evaluation rules (innermost to outermost) cannot be used. This command can be written in functional form as

`assuming`(F(G(e1, e2)), [a1, a2])

`assuming` is a procedure written in Maple; you can read its code with showstat. Although it's quite complicated code in general, you can see at the top that its first parameter is declared ::uneval. 

Returning to your original example, if for some reason you actually want to use the assumption only for the simplify and not for the odetest, and you want to do it in a single expression (i.e., without the intermediate assignment to res), use

simplify(odetest(mysol, ode) assuming []) assuming x>0

If you prefer, the [] could be replaced with nothing or anything or x::real or x::complex. The inner assumptions then supercede the outer when the inner expression is evaluated.

I think that the following will work in Maple 18. Some minor simplifications of the input syntax are available in more-recent Maple. My way of entering the system highlights that it is a matrix-vector linear system. There are other ways of entering this that others may consider simpler, although they do involve a bit more typing.

Is:= i__||(1..2);
sys:= {
    seq(diff([Is](t), t) =~ <1/2, -3; 2, -6>.<Is(t)> + <5*exp(-2*t), 0>),
    ([Is(0)] =~ [1,-1])[]
};
Sol:= dsolve(sys);

The exact solution that you get here looks more complicated than
it really is. That's because the coefficients are roots of
quadratic equations containing square roots. These numbers 
can be put in decimal form via:

evalf(%);

You also mentioned Maple Calculator. Many commands available in Maple are also available in Calculator. I suspect that dsolve is one of them, but I'm not sure. It's an extreme amount of work to enter alphabetic commands through the keyboard available in the Calculator. On the other hand, Calculator can be used to take a photo of the problem and translate it into Maple syntax which you can then upload to your Maple 18.

Name the two equations from (19) eq1 and eq2. Exactly as you've already entered them is fine; just begin the lines with eq1:= and eq2:= . There's no need to enter the next two equations; the solve command does that work:

Sol:= solve({eq1, eq2}, {U__mn, V__mn});

Then extract the coefficients like this:

ABs:= {
    seq(a[k] = coeff(eval(U__mn, Sol), W__mn, k), k= 0..2),
    seq(b[k] = coeff(eval(V__mn, Sol), W__mn, k), k= 0..2)
};

If you wish to assign the values of a[1], ..., b[3] (I'm not saying that it's a good idea to do this!!), do

assign(ABs);

When it's possible to solve a problem symbolically, it's usually useful to do that and then substitute the numeric values for the parameters. The symbolic solution provides invaluable insight into the physics of the problem. This problem can be solved symbolically by nearly any method that can be used numerically, including those posted by Kitonum and Tom. And here's another:

LinearAlgebra:-LinearSolve(
   Matrix(scan= band[1,1], [[280], rho*h*omega^2+~[100,444], [200]]),
   -C*<555, 6665>
);
                    [       /           2             \]      
                    [   5 C \111 h omega  rho - 217316/]      
                    [ - -------------------------------]      
                    [                 %1               ]      
                    [                                  ]      
                    [      /            2             \]      
                    [  5 C \1333 h omega  rho + 102220/]      
                    [- --------------------------------]      
                    [                 %1               ]      
                      2      4    2              2            
               %1 := h  omega  rho  + 544 h omega  rho - 11600

(A,B):= seq(eval(%, [h, C, rho, omega]=~ [3/100, 1, 7800, 222]));
                               -63994265     -549030931 
                     A, B := -------------, ------------
                             1330038150364  950027250260

evalf([%]);
                    [-0.00004811460858, -0.0005779107187]

Is there any indication in the documentation that algsubs might be appropriate for this?

Just use eval(ode, sol). Note that y(x) is necessarily a syntactic subunit^[*1] of any expression in which it occurs, so there's no reason to think that algsubs might work better; and, as always with algsubs, there are many  reasons to think that it might work worse.

Why do so many people gravitate towards algsubs? and then after algsubs, they try applyrule? That those commands are houses built on sand can be guessed just by reading between the lines of their help pages. And if that doesn't convince you, read their code. The fundamental commands are subs and eval.

[*1] If you understand the representation of expressions as trees or directed acyclic graphs (DAGs), the syntactic subunits are the nodes of those graphs.

Some of the below may duplicate things that acer already said. I didn't follow the links to worksheets that he gave (can't do it from my phone), so I haven't read the entirety of his Answer.

You wrote:

•  fsolve and plot only work with the proc if a function was defined in their calling sequence as shown below.

That's totally a red herring reinforced by confirmation bias. (Not your fault: The human mind evolved to fall prey to confirmation bias.) The truth is that in the examples that you show that didn't work, the first argument was not a proc; and in the examples that you show that did work, it wouldn't have mattered whether the proc was defined inside or outside of the command call.

Here are your examples (in 1D input), 1st those that didn't work. Given: f is exactly as you've defined above; it is a proc.

fsolve(f(x)[1] - x, w);
plot(f(x)[1] - x, 0. .. 0.02);

The problem is that f(x)[1] - x is not a proc. So the fact that f is a proc is made irrelevant in these examples.

Now those that worked:

fsolve(x-> ValveLaminarMassFlowProc(x, P__1, P__2, rho__1, T__1, rho__eta)[1] - x, w);
and likewise for the plot.

You could just as well have defined:

f1:= x-> 
    ValveLaminarMassFlowProc(x, P__1, P__2, rho__1, T__1, rho__eta)[1] - x
:

and then done: 

fsolve(f1, w);
plot(f1, 0. .. 0.02);

and those would've worked also, which proves that having the procedure definition inside the command is irrelevant.

You also asked about trapping the error inside the proc. I suppose that that can be done, but I'd rather show you how to prevent the error inside the proc. This will involve two programming concepts that may be new to you: type-checking and returning unevaluated.

Given: The same f and f1 as above, define:

f2:= proc(x::algebraic)
    if x::numeric then f1(x)[1] - x else 'procname'(x) fi
end proc:

The x::numeric is the "type check", and the 'procname'(x) being the last thing evaluated when the else branch is taken is the "returning unevaluated". The keyword procname is used (only inside procedures) to refer to the proc's own name (f2 in this case), i.e., it's a form of indirect self reference. The quotes around the procname are what makes it unevaluated.

That's all there is to it! It's not actually necessary to have a cascade of procedures f, f1, f2. I just did that to save myself some typing. All this can be done in the original f.

Now you can use f2 like this:

plot(f2(x), x= 0. .. 0.02);
and likewise for fsolve.