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.

My own recent bad restart experience:

What you describe is a serious bug in restart. I also had a recent (about 3 weeks ago) experience that convinced me beyond any reasonable doubt that

1. Some kernel memory had survived a restart,
2. That memory was causing an error, and
3. An expression that had been part of the computation before the restart was being blamed in the error message as the cause of the error even though the computation after the restart did not have the remotest connection to the one before.

When run in a new kernel, the second computation ran without error (as expected). I was able to duplicate this whole scenario multiple times. The expression was a totally ordinary algebraic expression that had been the denominator of an ODE IVP that I had solved numerically (via method= classical[rk4]). Indeed, the expression was d:= (x(t)^2 + y(t)^2)^(3/2). It's unique enough that there's no chance that it would appear anywhere at all in the second computation simply by coincidence. The second computation had nothing to do with differential equations or even algebra. The command that errored was (I think) type(e, with_unit(realcons, 'a', 'b')) where e was something very simple like 100*Unit(ohm). (I'm sure about the unit, but not about the coefficient. The coefficient was definitely some positive real number.)

Minutiae regarding my enviroment (in case it may help someone investigating this):

1. Neither computation involved any user-initiated disk I/O (e.g., save, read, savelib) whatsoever, nor were any directory-checking or directory-changing commands used (e.g., currentdir).
2. I almost never load packages with with, and I certainly didn't in this situation.
3. I use Standard interface, 1D Input, 2D output with extended typesetting.
4. Windows 10; Maple 2020.1 (I don't trust .2).

5. My initialization file is very basic. My libname includes a few standard 3rd-party packages such as DirectSearch, but none of those were used in either computation.

6. I don't get the Physics updates because installing them seems like a hassle and I do no symbolic work in physics and very little in special functions or differential equations. (I mostly do numerical analysis, combinatorics, operations research, symbolic linear algebra, logic, and stuff with very complicated types.) However, this new Latex is intriguing, and I hope that it's preinstalled in the next major release.

7. I used no embedded components in these two computations, not even a plot.

(I have no idea why the MaplePrimes editor double-spaced those last two points, which I've tried to correct several times.)

What restart is supposed to do:

First: (You probably know this already, because I've said it many times over the years (not necessarily to you).) You need to put restart in its own execution group. This is mentioned on its help page. However, it does usually---but not always---work correctly if you don't follow this advice.

AFAIK, restart is supposed to do all these things:

1. Execute the ModuleUnload procedure of any module in use that has one. (Most modules do not, but I think that Physics has a very long one. You can see what it's doing by setting printlevel:= 10 (or greater) in a separate execution group immediately before doing the restart.)
2. Reset all variables to their initial state.
3. Collect the garbage.
4. Execute your initialization file, if you have one.

And I'm unsure about whether it's supposed to do these:

1. Close files opened by the user. Surely, it should close them (after posting their write buffers), because their pointers will disappear with the restart, thus making them zombies. But I don't see documentation that it does close them.
2. Update the status bar to reflect the fact that all the garbage has been collected. (I'm sure that it doesn't do that updating; the part that I'm unsure about is whether it's supposed to. It seems shoddy that it doesn't.)
3. Things that it's supposed to do, does do, or doesn't do to the GUI are of no concern to me at the moment (except for that status bar point).
    

Something that you can do as a workaround that's a lot faster than what you're doing now:

If you save and reopen your worksheet, that will start a brand-new kernel (a new invocation of mserver.exe). This takes much less time than opening and closing your entire Maple session.

And, as I've said several times before (not necessarily to you), this exact same advice can be used to correct loss-of-kernel errors (rather than closing the entire Maple session, as is incorrectly recommended by the loss-of-kernel error message).

Yes, it's simply indices(node).

P.S.: I believe that an earlier version of this Question referred to node as a "Matrix". I guess that you realized that a matrix can only have two indices and thus edited out that word. So, let me reassure everyone that the command indices will work exactly as shown on any table or rtable (Array, Matrix, or Vector).

You are ignoring everything that I tried to teach you in my responses to your previous version of this Question, thus causing unnecessarily duplicated effort for other responders. Here's a summary:

1. Separate the code that numerically evaluates the function from the code that numerically integrates it.
2. Replace your conditional (if-based) definition of r with an unconditional one using abs. You already agreed, after some cajoling, that this is correct, so why go back to the if?
3. There is some severe numerical instability (or perhaps it's extreme oscillation) in this function that has thwarted every attempt that I've made (using up to 8 million evaluation points) to numerically integrate it (using many different algorithms) to more than 1 digit precision. Your attempts to use the trapezoid rule or Simpson's rule (here and in the prior Question) are worthless if you can't get an upper bound on the error.

Do you think that all the advice that I already gave is worthless just because it didn't finally produce the integral?

Do you think that adding up a bunch of nearly random numbers and calling it an "integral" is worthwhile? 

To proceed in a meaningful way, the numerical error needs to be controlled. I'm hoping that some other respondent can provide some insight into that. Once the error is controlled, one of Maple's standard numerical integration algorithms (which are far more sophisticated than Newton-Cotes (trapezoid, Simpson's, etc.)) should be able to handle it. That makes my point #1 above even more relevant.