sys := {z = 4, x-y = 5, x+y = 10};
L := map( rhs, solve(sys, {x, y, z}));

What is wrong with

fq := p -> piecewise( p < -5, 0.,
                                   p < -3, .161711971*(p+5.)^3,
                                   p < -1, -1.798464149*p^2-.544033315*p+7.541565235-.3076373305*p^3,
                                   p < 1, 0.809338641e-1*p^3-.6327505649*p^2+.6216802691*p+7.930136430,
                                   p < 3, 0.7445052690e-1*p^3-.613300552*p^2+.602230257*p+7.936619768,
                                   p < 5, 0.19182821e-2*p^3+0.39489651e-1*p^2-1.35614035*p+9.894990378,
                                   p < 7, -0.256603524e-2*p^3+.106754411*p^2-1.69246414*p+10.45553002,
                                   p < 9, -0.254356370e-2*p^3+.106282509*p^2-1.68916084*p+10.44782231,
                                   p < 11, -0.2031067084e-1*p^3+.585994401*p^2-6.00656786*p+23.40004343,
                                   p < 13, 0.4137545813e-1*p^3-1.449647855*p^2+16.38549696*p-58.70419431,
                                   p < 15, -0.273325019e-1*(p-15.)^3,
                                  15 <= p, 0.
                                ):
m := p -> piecewise( p < 1, undefined,
                                   p < 9, fq(p),
                                   undefined
                                ):
plot(fq, -20..20);
plot(m,  -20..20);
seq(fq(j), j=-20..20);
seq(m(j),  j=-20..20);

If OP just wants a simple verification then

restart;
with(Logic):
Implies( (a &implies b) &and (b &implies c),
                a &implies c
            );

will provide it. (May be redundant parentheses in the above, logical operator precedence is not my strong suit, so I preferred to be safe rather than sorry)

If you can guarantee that the only arguments you want are strings, then

restart;
propositionsentence := MM(AA("it", "run"), "maple"):
indets(propositionsentence, string);

will work - it returns {"it", "maple", "run"}.

If your test case contains a mixture of strings and names and you want both, as in

propositionsentence := MM(AA("it", "run",x^2), "maple", sin(y)):

then

indets(propositionsentence, {string, name});

will return

{"it", "maple", "run", x, y}

Leaving aside the fact that the posted command has the 'restart' command in the wrong place :

I can see two ways when you *might* get an odd number of roots

1. supplied value of h_given coincides with a local/global maximum of the function 'f'. You will be "unlucky" to hit this condition accidentally - but you might. Note that if you hit a local maximum, you will generate an odd number of roots, but the number of roots will be >1. In order to get exactly one root, the value of h_given would have to coincide with the global maxmimum of f
2. You search for roots in the range 0..timeT+5 (although you only plot the range 0..timeT). If the value of h_given happens to match the value of f in the range timeT..timeT+5 - then it is possible to produce a single root. So again you will generate an odd number of roots. In order to generate just one root this way, the function f would have to attain a value between timeT and timeT+5, which it does attain between 0 and timeT. Very unlikely, but I guess it is theoretically possible. You can remove this possibility by change the search range in the Roots() command to 0..timeT

I can produce an answer, but I'm not sure if it is what you want/expect.

In the attached

1. I defined the function epsilon(rho) as a piecewise function over the three ranges you requested (and undefined outside these ranges)
2. In one of your odes you use epsilon, rather than epsilon(rho), so I fixed this
3. You supply four BCs, but you can only use three - so I just (arbitrarily) threw one away
4. Maple will now treat the problem as a BVP over the range r=0....1.0e-05, so one can only get solutions within this range
5. In order to obtain a numeric solution, one has to supply a value for 'rho' in order that epsilon(rho) can be calculated. One could use multiple values for 'rho' in some kind of loop - depends on what format of resultss you want. I don't know what you are trying to achieve here so I haven't bothered (and in any case results don't seem to vary very much as 'rho' is varied, so I haven't bothered)

Check out the attached to see if it helps

odeProb.mw

fsolve(diff(fse32_2, p) = 0, p = 1 .. 9)

finds the value of 'p' at which the maximium occurs so fast, I couldn't be bothered to time it!

the PDETools:-Infinitesimals() command??

It is supposed to worrk for either a single PDE or a system of PDEs

I used the add-on package DirectSearch which returned the solution

c[0][0] = 12.8960498385722,
c[0][1] = -6.47066667098714,
c[0][2] = 1.75211104639280,
c[0][3] = -.177494271892552,
c[1][0] = 9.73759571914550,
c[1][1] = -2.59015291121796,
c[1][2] = .340246167875447,
c[1][3] = -0.170026496357586e-1,
c[2][0] = 7.03216883026407,
c[2][1] = -1.17901185303931,
c[2][2] = 0.969676292786482e-1,
c[2][3] = -0.312131696475318e-2

This agrees with second (real) solution you obtained with solve() and the solution Preben obtained with fsolve(), when provided wth a reasonable guess

Despite playuing around with a few options, (so far) I haven't been able to persuade DirectSearch to produce the first (real) solution which you obtained from Maple's solve() command. Not saying it can't be done - just that I haven't been able to do it :-(

You seem to hav a fundamental misunderstanding about what a contour plot actually is!

Now Kitonum's earlier answer is correct, but I am going to try to break things down into simpler pieces for you.

First of all, I assume that you have no difficulty with

f := exp(-r/2)*r^2*cos(2*phi);
plot3d(f, phi = 0..2*Pi, r = 0..15,  coords = cylindrical);

Now I want you to look at this plot - very carefully. In particular, when (in cylindrical coordinates) z=15, you have a four-lobe structure, which extends over ( very roughly, by eyeball) about x,y=+/-0.2.

So try plotting the single contour corresponding to f=15, using

p1:= plots:-contourplot( f,
                                          phi = 0..2*Pi,
                                          r = 0..15,
                                          coords = cylindrical,
                                          contours=[15],
                                          color=red
                                       );

Now this will clearly show that in fact the lobes on the original 3-D plot only extend to x,y=+/-0.15 (approximately)

Similarly the 'fattest' point on the original plot3d graph is at z=4 (roughly), where x,y values extend to +/-2 (roughly), which you can see with

p2:=plots:-contourplot( f,
                                         phi = 0..2*Pi,
                                         r = 0..15,
                                         coords = cylindrical,
                                         contours=[15],
                                         color=blue
                                      );

Now it is possible to put both of the above contours on the same graph, using

plots:-display([p1,p2]);

However if you use the command

plots:-contourplot( f,
                                 phi = 0..2*Pi,
                                 r = 0..15,
                                 coords = cylindrical
                               );

then Maple will provide (by default) eight contour curves, equally spaced between your minimum and maximum z-values, ie 0 to 15, so the spacing will be 2143 (approximately)

On the other hand you might want to control the spacing of the contour values. An obvious choice for your problem woul be to plot the contours for z =0,1,2,3,,,,15, which would be obtained with

plots:-contourplot( f,
                                 phi = 0..2*Pi,
                                 r = 0..15,
                                 coords = cylindrical,
                                 contours= [ seq( j,
                                                             j=0..15
                                                           )
                                                   ]
                                );

HAve a long hard look at the attached

contourStuff.mw

You can use the 'events' option on the dsolve/numeric command, which will force some action to be taken when a specified criterion is met during the dsolve process.

In your case the criterion is ca(t)-0.2=0 - and the simplest 'action' to take when this occurs is just to halt the dsolve() process, which will cause the values ca(t) and 't' to be output when ca(t)-0.2=0. As in

restart;
sys:= diff(ca(t), t) = -3.600000000*10^20*exp(-15098.13790/(340-20*ca(t)))*ca(t), ca(0) = 2:
sol:= dsolve({sys}, numeric, events=[[ca(t)-0.2=0, halt]]):
interface(warnlevel=0):
sol(10);
interface(warnlevel=3):

which returns

[t = 3.89994725315188, ca(t) = .200000000000000]

In the above, the sol|(10) command will attempt to compute the solution for 't=10'. In order to do this, it has to start at 't=0', and then iterate its way to 't=10'. However during this process, if it ever sees the condition 'ca(t)-0.2=0', it will 'halt' and output relevant values at that point.

By the way, you don't really need to toggle the 'warnlevel' parameter as I have done in the above - Maple will issue a warnng whenever the dsolve() command is terminated by a requested 'event' condition. Since this 'event' is what you/we are looking for, it seems (to me) odd that Maple throws a warning. So I have a tendency to turn  'warnings' off when I am doing something like this, and then turn them back on again (at the default level) when I have finished

As noted by Kitonum, if you plot all curves on the same graph, you will only 'see' the final section of the quadratic plot - everything else will be "collapsed" on to the x-axis.

So I plotted the three curves separately

The following Kitonum's suggestion I used plots:-logplot, to make all three curves vsible on the same graph

someplots.mw

The attached illustrates the immediate answer to your problem, together with a few (non-onbvious) ways in which things may go wrong. I would consider some the latter to be poorly or completely (un)documented.

procArgs.mw

Make of it what you will

You cannot split the do loop across multiple execution groups.

At the Maple prompt, start typing your for loop, and use Shift+Enter for each (soft) new line. Keep typing till you have entered the 'end do' statement. Then use "Enter"

So weird, in fact, that I ended up more or less rewriting your code in my attempt to figure out what is going wrong!

In the attached,

1. I generated your specified 512 matrices, and put the corresponding digraphs in a list
2. For each of the digraphs in the list, I check whether it is isomorphic to every other digraph in the list
3. I use a' try-catch' construct to ensure that execution doesn't just halt, when an error occurs, and I can save the list indices for which the isomorphism test fails
4. Failure is always associated with error message "Error, (in GraphTheory:-IsIsomorphic) invalid subscript selector", and I always get three failures, associated with the same indices into the list of digraphs
5. I can delete the list of digraphs, regenerate it from scratch, and check only those indices which failed in the global check. They still fail!!!

See the attached, where I have have saved the worksheet with output so you can see the digraph indices for which failure occurs.

weird.mw

This error is so weird that it might be OS or Maple version specific, so if you rexecute this code, you might get different answers.

The error message is incorrect, because there is nothing wrong with any subscript over which I have control. I am 99% sure that this is a Maple *bug*, and I can't figure out how to work around it :-(