@Carl Love After reading the example shown by dharr, I realized that each of the eigenvectors corresponding to eigenvalue 0 has with itself precisely the type of symmetry that you want. These eigenvectors are always of the form < x, -x >.  Let < u, v > = < x, -x >. Any scalar multiple of an eigenvector is also an eigenvector for the same eigenvalue, so - < x, -x > = < -x, x > = < v, u >, corresponding to eigenvalue -0 = 0.

This is fairly easy to prove rigorously. (I'd consider giving it as an exercise to an honors undergraduate linear algebra class.) Let me know if you need help with that.

So, the symmetry anomaly that you're seeing is not related to the numerical anomaly. And I think that the numerical anomaly will resolve when you specify shape= symmetric.

@jefryyhalim Yes, I can reproduce the erratic behavior as you described. My first reaction is to increase Digits to 30. I start getting undefineds on line sol_L:=.... By the maxim I just stated, the results for all lower values of Digits are unreliable.

See how far you can get doing exact computation only. Convert all decimal numbers in the input into their equivalent fractions. Now, this may quickly get bogged down with expression swell. Or it may work. I recommend that you execute line by line. Remove all colon statement terminators so that you can observe the output of every line. As soon as you see a decimal anywhere in the output, correct the corresponding input to exact values.

I am suggesting this, at the moment, merely as a means to find the error (if any) in your logic. I'm not suggesting, at the moment, that the final production code will need to be done in this exact form.

@jefryyhalim Here's a maxim for all numeric computation (not just Maple): Never ignore the results provided by using a higher value of Digits when they differ strangely or significantly from those provided by a lower value unless an expert in numeric computation has looked over your problem and given you a specific computational/mathematical reason to ignore them. In the vast, vast majority of real-world cases (I'd rough guess at least 99.9%), the results from the higher value of Digits are closer to the truth. It may be a truth that you don't want to see, but it shows that the results from the lower value of Digits are completely unreliable. To ignore that is putting your head in the sand. That it makes the model fit the experimental data is never a sufficient reason. 

@mmcdara What the OP has written is correct when used in 2D Input. To verify this, copy the OP's input line V(r):= ... to a 1D Input (aka Maple Input) execution group. You'll get this (in Maple 2018):

V := proc (r) options operator, arrow, function_assign; piecewise(0 <= r and r <= R, -VV, R < r, 0) end proc

I think that the simulation that you want can be done in Maple itself. If you post details and formulas, we can get started.

@jefryyhalim I wouldn't say at this point that the behavior arises from the ODE. The red flag to me at the moment is this 1e-10, which I'll call epsilon.

To be clear about what we've learned so far: You know now that you can't use Digits = 10 with this small a value of epsilon, right? Digits:= 15 is a very good value to use because it's the largest value such that computations are done in the realm of hardware floats, which are numbers for which the very fastest computations are possible.

Now, you need to explain to me What is your purpose in using epsilon? I find it quite suspicious, and a likely cause of erratic, unstable, or undesirable behavior. I notice that you divide things by tan(Pi/2 +/- epsilon), which is essentially +/-infinity. This is exactly the kind of thing that is likely to cause Float(undefined) to appear.

@jefryyhalim Early in the worksheet, you set a variable to (essentially) evalf(Pi/2 + 1e-10). If Digits is 10 (its default value), then this is exactly equal to evalf(Pi/2). Surely that wasn't what you intended. I see that you had Digits:= 15, but you commented it out. Why? If you use Digits:= 15, do you still have the erratic behavior?

@max125

Yes, the purpose of the subs(y(x)=y, ...) is so that y(x) will be treated as a regular variable (aka, a name).

Yes, without any knowledge of the nature of the function y(x), it is impossible to solve any equation containing y(x) for x. It may still be solved for y(x) or another variable.

The process of treating a function (or other expression) as a simple variable is common enough that there's a command just for it: frontend. In this case,

frontend(solve, [deq= 3, y(x)])

works. Its default action is to treat all functions appearing in the square brackets as simple variables. In this case, the only function is y(x). This can be controlled with options. See ?frontend.

Meaning of RootOf: If f(_Z) is any algebraic or transcendental expression containing variable _Z (and possibly other variables), then RootOf(f(_Z)) means "the set of all _Z in field F such that f(_Z) = 0." Here, F represents the algebraic field of interest. Usually this is the complex numbers, although Maple can deal with other fields as well. Maple often responds with a RootOf  when it can't explicitly determine the members of that set. There's only ever one bound variable in a RootOf, and it's always _Z.

Start each repetition of the computations with the command

restart;

by itself, on its own line, in its own execution group. Beginning with restart is good advice for any Maple session. If you still experience the erratic behavior that you described, please repost here including an upload of the worksheet.

On Github, there's a probabilistic programming language called Hakaru that's written partly in Maple and partly in Haskell.

What do you mean by "add support for" non-plaintext worksheet files .mw and .mws? I understand what "support" for plaintext files like .mpl means. And if you're going to do worksheet files, you might as well also include .maple, the newest member of the family.

@rahinui There are 4 types of remember tables that can occur in Maple procedures (if you count absence of a table as a type). This is mostly^[1] controlled by the options (if any) at the beginning of the procedure. Every procedure has exactly one of these types:

• Regular (option remember): All computed values are remembered until they are explicitly removed.
• System (option remember, system): Just like the above, but the table is also cleared by gc() (whether automatic or explicit).
• Cache (option cache): The table has a maximum number of entries (that you can set, if you want). When it's full, new entries replace the oldest ones. These tables can also have "Permanent" entries, which are explicitly set.
• No remember table (there is no corresponding option): The procedure is executed whenever it is called.

I am considering a cache table a type of remember table. The documentation seems to distinguish these, but in name only: There's not much difference between them (other than as I described above) as far as I'm concerned.

Regardless of the options, a procedure doesn't acquire a remember table until there's actual data to store in it.

Regardless of the options (or their absence), entries can be added to the table of procedure P using the syntax P(...):= ..., where the "..."s represent any expression or expression sequence, including NULL. If the procedure doesn't have a table, one will be created. If the procedure doesn't exist, this will create it.

This procedure checks whether a procedure has a remember table:

HasTable:= (P::procedure)-> evalb(op(4, eval(P)) <> ()):

And this one checks whether a procedure has the option for one:

HasTableOption:= (P::procedure)->
   evalb(({op(3, eval(P))} intersect {'remember', 'cache'}) <> {})
:

^[1]"Mostly" because there's also direct-assignment: P(...):= ....

(to be continued)

@Arny

Surely you're closer to the result than you were before I intervened!?

I may not being understanding the "Physics" aspect of this correctly, but how could the result possibly be anything other than infinity? The integrand doesn't even approach 0. Simplify it to the 1D case, and perhaps you'll see that. I'm suggesting this as a "thought experiment" at the moment; you still have some syntax issues that we can address later.

@Britzel Could you give me a simple example of a functional F of a single argument f (itself being a function of a single real variable x on a closed interval) that you might want to minimize?

@Christian Wolinski First vote up is mine.

Just to clarify for the other readers: What you're using here is called simplify with side relations, and it's documented at ?simplify,siderels. The factor and the variable-order argument seem superfluous (in this particular case at least), and simplify automatically maps itself, so your code could be

simplify(
   evecA1[2],
   [rho*c^2/(gamma-1)*(1/rho)+(1/2)*u^2 = H]
);

The command algsubs is much overused. I find its results less predictable than its more-modern alternatives such as the above and subsindets (which won't work in this case).

@Mariusz Iwaniuk I don't think that you understand Joe's comment about the op in your Answer. For the particular example set or list given in the Question, the op(1 , ...) does nothing, and thus your code gives the correct result. But if a minor change were made to the example, such as  L:= {a = b[1], e = f[2]}, then that code would give an incorrect result.