Set 

kernelopts(numcpus= 1)

Many years ago, I discovered a bug such that one could only change kernelopts(numcpus) immediately after a restart. If it was changed afterwards (I guess after the first garbage collection occurs), the kernelopts command would claim that it was changed, but in reality it wasn't. I don't know if this has been fixed.

You wrote:

• My code does not do any mutlithreading. So I do not need it.

You don't know that! You do very high-level stuff that calls thousands of library procedures. Some of those probably attempt multithreading if they see CPUs available for it. But, I've never seen any Maple library command error out due to the absence of multithreading. 

Here's a procedure for it. It returns a domination set of minimal size so that its dominance can be easily verified. Then it's trivial to count it if you want.

Download DominationSet.mw

You wrote:

• Since I cannot feed non-algebraic expressions into polynomial solvers, I want to extract only their numerators (which should be algebraic)....

You are misusing the mathematical term algebraic. Since I think that there's a good chance that such misuse will lead you to make even-more-serious errors, I'm going to digress for a bit to attempt to correct that misuse. First, note that Maple has a type named algebraic. That type has nothing to do with the mathematical concept! This discussion is strictly about the mathematical concept.

I'll assume that you understand the definitions of all these mathematical objects; let me know if you need clarification on any:

• integer
• rational number
• polynomial
• polynomial in the variables {x, y, etc.}
• polynomial in the variables {x, y, etc.} with integer coefficients
• polynomial in the variables {x, y, etc.} with coefficients that are polynomials in the variables {a, b, etc.} with integer coefficients.

Your terminology usage suggests that the above is perhaps the limit of your knowledge of the precise definitions of even-larger classes of objects. Thus:

Definition: A rational function is a quotient of polynomials, or an expression that can be re-expressed as such via reversible common arithmetic operations (as learned in high-school algebra). As above, we may refine the class by specifying the variables and the type of coefficients.

Definition:  An algebraic function is a rational function with the possible inclusion of non-integer rational-number (fractional) exponents (such as square roots). These exponents may appear not just on individual variables but also on any subexpressions.

Let me know if you need any clarification of those definitions.

Now, returning to your immediate questions, it seems that by "non-algebraic" you actually mean non-polynomial.

• These solutions should also solve the non-algebraic [non-polynomial] system then.

Almost true. The solutions thus obtained may be a nontrivial superset of the solutions of the original system. They need to be verified, which is a much much easier process than obtaining that superset. In the case at hand, your original system is rational functions, and thus the only thing that really needs to be verified is that the solutions do not make any of the original denominators zero.

• To extract the numerators, I am using ((numer@evala@:-Norm@numer)~@eval)(Eqs)....

Obviously, that's either a command that I composed, or a variation of such. In this case, all that you need is numer~(Eqs) because they are rational functions. And the eval is only needed if you want to instantiate parameters. But deconstructing anyway:

1. @ is the function-composition operator: (f@g)(x,y) = f(g(x,y)). As you may recall from a long-ago math course, function composition is associative but commutative. Let me know if you need a refresher.
2. ~ is the elementwise operator (see help page ?elementwise). It causes the operator to which it is appended to be applied individually to the elements (members) of a container (such as a set) that follows. In this case specifically, it applies the operators to your individual equations. (By equation of course I mean an expression implicitly equated to 0.)
3. evala is a "package" (in an older Maple sense of "package") of commands for performing sophisticated operations on algebraic functions. Since you only have rational functions in this case, it's not needed. When I originally gave you the command, the situation was much more complicated with numerous radicals in the equations. Those were algebraic, but not rational, functions.
4. Norm is a command in evala for converting algebraic functions to rational functions. This process necessarily involves irreversible operations (such as squaring both sides of an equation to remove square roots). Of course, it's not needed in this case.
5. :- is a purely syntactic "disambiguation" operator. On the off-chance that you had loaded some other package (such as LinearAlgebra) that also has a command named Norm (totally unrelated), the :- specifies that that is not the Norm that we mean.
    

• 1. Am I wrong?

No. Other than the terminology issue and using a needlessly complicated command for the easy job of extracting the numerators of rational functions, you're totally correct.

• 3. Why length(((numer@evala@:-Norm@numer)~@eval) numer~(Eqs)) is much larger than length(Eqs)? Shouldn't it be smaller since I just extract the numerators?

Example: Let e:= 1/x + 1/y + 1/z. The length of e is 25. The numerator of e is y*z + x*z + x*y. It's length is 37. Although it's not important at this point to understand how the exact values 25 and 37 are obtained (it's idiosyncratic to Maple), isn't it obvious that the numerator alone should be considered "lengthier" than the original? If it isn't obvious, extend e to more terms and see what happens.

Finally, your title also mentions "simplify". Other than instantiating parameters to numeric values, I wouldn't simplify, unless absolutely necessary due to memory usage. Compared to the solving process, it's trivial for the polynomial solver to simplify in whatever way works best for it.

If you don't instantiate at least some parameters, I think there's very little hope of obtaining any solutions. I'm not saying that those solutions don't exist, just that I think that they're not possible to obtain symbolically.

Your procedure F does not make any assignment. In particular, it makes no assignment to f. Thus, no matter how many times you call F(i), and no matter what value you use for i, the only assigned values of f will be f[2] and f[3]. To correct this, the statement in the for loop should be

    f[i]:= F(i)

If you read the extremely simple and extremely short code of MmaTranslator:-Mma:-LeafCount, then I think that you'll immediately see what's going wrong. Basically, it attempts to break everything down into numbers and names. Your input is a string, so it cannot be broken down into numbers and names. If you remove the quotes from your input, then you'll get the 1 that you expect. If you still can't understand it after reading the code, let me know, and I'll explain it in more detail.

The entire Mma subpackage of MmaTranslator reads to me as if it was written by a rank amateur, such as someone who had 1 or 2 programming courses freshman year at university, then maybe self-studied Maple for a month.

There can only be one "plot" per file. I put "plot" in quotes because any number of plots can be combined into a single plot by passing an Array of plots (rather than a sequence of plots) to plots:-display. This can be done simply by enclosing the seq command in angle brackets: <seq(...)>.

Example:

plotsetup('png', 'plotoutput'= "ArrayOfPlots", 'plotoptions'= "width=1920,height=1080"):
cols1:= ["red", "green", "blue"]:
plots:-display(<seq(plot([$20], 'color'= cols1[j]), j= 1..3)>);
plotsetup('default'); #Reset output to normal screen

I don't know whether the above will be useful to you. If you need to afterwards separate the individual plots, then it may not be very useful. In that case, you can create essentially an "array of files", 1 plot per file, like this:

restart:
cols1:= ["red", "green", "blue"]:
Plots:= <seq(plot([$20], 'color'= cols1[j]), j= 1..3)>:
PlotToFile:= proc(P::specfunc(PLOT), j::posint)
    plotsetup(
        'png', 
        'plotoutput'= cat("ArrayOfPlots_", j), 
        'plotoptions'= "width=1920,height=1080"
    );
    (* Threads:-Sleep(2); *) #Sometimes helps
    print(P)
end proc
:   
seq(PlotToFile(Plots[j], j), j= 1..numelems(Plots));
plotsetup('default'); #Reset output to normal screen

I didn't need to use the Threads:-Sleep. Sometimes it's needed to let a slow disk drive catch up with the faster output from Maple. The 2 means to pause 2 seconds.

I appreciate and understand what you're advocating regarding a programming paradigm that avoids environment variables. In addition to what you said, they also hinder attempts at shared-memory parallel processing (as offered by the Threads package) when different values are used in different threads. Unfortunately, their use (as well as non-environment global variables) is endemic in the Maple library code, which is I guess the main reason why there are so few commands listed at ?threadsafe. Just to remove the usage of the single environment variable Digits would be a monumental undertaking I think. I'd guess, very roughly, that about a thousand library procedures use it.

I think that there are only two differences between environment variables and other globals:

1. Environment variables don't need to be declared in procedures.
2. When an environment variable's value is changed in a procedure, it reverts to its previous value when the procedure exits.

Any global variable whose name begins _Env can be used as an environment variable. And there are a few very common ones such as Digits and Order whose names don't follow any special pattern.

_EnvTry doesn't have an initially assigned value. That's why you don't see it listed by anames(...). Your procedure is correct. But note that when the procedure ends, it'll revert to not having an assigned value. That's not a problem. I just mention it so that it doesn't come as a surprise to you.

1. Rotation: The wobbly rotation that you want can be easily done with plottools:-rotate applied to the static plot3d p that you already had. The [[.1,.1,-1], [.1,.1,1]] is the line segment parallel to the z-axis to rotate about.

plots:-animate(plottools:-rotate, [p, B*Pi/36, [[.1,.1,-1], [.1,.1,1]]], B= 1..72, frames= 72)

2. Field lines: I think they'd be better done with plots:-tubeplot because that makes them surfaces thus giving you a much larger set of options (shading, color, glossiness, transparency, etc.).

3. Aesthetics: axes= none. Set thickness=0.2 to make the gridlines (on sphere) thinner. But if you remove them entirely, it makes the rotation more difficult to see IMO.

If you plug the solution (converted from series to polynom) back into the ODE and simplify:

simplify(eval(lhs(ode), convert(sol, polynom)))

then all that remains are terms with orders equal to or greater than the O(...) terms in sol. That's a reasonable indication that those "remainder" terms are approaching 0 as Order approaches infinity. 

Maple considers a RealRange to be a property rather than a set (in the Maple sense of set). The equivalent of union for properties is OrProp (always used in prefix functional form: OrProp(RealRange(..., ...), ...). The boolean evaluator for this is is, for example,

is(1/2, OrProp(RealRange(0,1), RealRange(2,3)));
                              true

Assuming that you want to start the indices at 0:

add(add(a[i, j, 2 - i - j]*x^i*y^j*z^(2 - i - j), j= 0..2 - i), i= 0..2)

The reason why solve acts differently in this case: Since isolate only makes sense for a single equation, the square brackets on the first argument are never necessary; so they weren't included in the allowed syntax. On the other hand, solve can solve multiple equations, so square brackets are allowed.

@mmcdara You asked:

• A question for more skilled Maple users: Is it true that dsolve/numeric is not capable to handle a Fourier boundary condition for a first order ode?

I don't know what "Fourier boundary condition" means, but I do understand what is happening here, and it's almost always easy to correct. It's not specific to first-order ODEs.

The issue is that dsolve(..., numeric) won't accept boundary conditions that use a derivative of the same order or higher order than the highest order of derivative of the corresponding function occuring in the ODEs.

Your example is

{diff(Phi(Y), Y)+f-Phi(Y) = 0, Phi(1)=1.5*D(Phi)(1)}

The boundary condition uses the first derivative of Phi, which is the highest-order derivative of Phi used in the ODE. Algebraically solving the ODE for that derivative yields diff(Phi(Y), Y) = Phi(Y) - f. So, in the BC, replace D(Phi)(1) with eval(Phi(Y) - f, Y= 1). The same thing can easily be done with the OP's system.

Here's a procedure for it:

RelErr:= proc(f::algebraic, V::list(name= `&+-`(realcons, realcons)))
local x;
    eval(add(abs(diff(f,lhs(x))*op(2, rhs(x))), x= V)/abs(f), lhs~(V)=~ op~(1, rhs~(V)))
end proc
:  
RelErr(R*i^2*t, [R= 1.5 &+- 0.5, i= 2.5 &+- 0.05, t= 10 &+- 0.5]);
                          0.4233333333

RelErr(G*M*m/r^2, [M= 7 &+- 0.05, m= 2 &+- 0.02, r= 5 &+- 0.05]);
                         0.03714285714

Use fsolve, not solve. Nothing else needs to be changed. The problem is that solve is trying to find all complex solutions.