There is no programmatic reason why you would need to add parentheses. If you want to add parentheses to expr for display purposes, do expr_par:= ``(expr). To remove those parentheses, which you'd need to do to use expr programmatically, do expand(expr_par).

Scot has given the answer. I just wanted to add some details. Items in { } are called sets in Maple. Sets are always sorted by Maple into an order that it finds useful (for efficiency reasons) that you have very little (usually none) control over. Items in [ ] are called lists in Maple. Lists stay in the order that you made them.

Like this:

num:= [a,b,c]:
den:= [c,b,a,d]:
sys1:= DynamicSystems:-PrintSystem(num, den);

In the above, Maple considers num and den to be lists rather than vectors, but there is no mathematical relevance to this subtle distinction. If for some reason you need them to be actual Maple vectors, you can simply change the above to this:

num:= <a,b,c>:
den:= <c,b,a,d>:
sys1:= DynamicSystems:-PrintSystem([seq](num), [seq](den));

The next thing to try is a continuation parameter. Read about it on help pages ?dsolve,numeric,bvp and
?dsolve,numeric_bvp,advanced. This is a parameter C that you include as a coefficient in your system. It will range gradually from 0 to 1. You place it such that when C=0 the system is easier to solve and when C=1 you have your original system. Note that the same parameter can appear any number of times in your system.

Like this:

J:= {$1..3}:  S:= {$1..2}:
Constraints:= {
    seq(seq(add(pi[i,j,r,F], i= J) = 1, j= J), r= S),
    seq(seq(seq(add(pi[i,j,r,u], i= J) = 1, j= J), r= S), u= S),
    seq(add(alpha[j,u], u= S) = 1, j= J),
    seq(seq(gamma[j,u] + add(gamma[j,r,u], r= S) = 1, j= J), u= S)
}:
simplify(Expression, Constraints);

I don't know what you mean by "all of the above-mentioned terms can be additively or multiplicatively connected"; I'll need to see your example of that.

Maple's default definition for a^b is one that makes the most sense over the largest regions of complex numbers possible. It's called the principal value, and an example of it is shown in nm's Answer. Under this definition, for a < 0 and b noninteger and real, a^b is never real. However, Maple allows for the redefinition (or overloading) of operators such as ^, and such overloading is often contained in packages. A package that overloads ^ to do what you want is RealDomain. So this will give you the plot that you want:

use RealDomain in plot(x^(1/3), x= -2..2) end use;

There's no specific limit, but large worksheets are a major problem. They cause the GUI to act very sluggishly. When the GUI becomes unresponsive, you must kill your entire Maple session (as opposed to killing a single kernel when a computation becomes unresponsive). Why not divide it into multiple worksheets? If the issue is that the worksheets must share data, that can be done. One way to share data is by collecting worksheets and data into a workbook.

Like this:

MyEqn:= (theta,k)->
local x, t:= theta*Pi/180;
    fsolve(x*cos(t) - (2500-x^3*sin(t))^0.45 - k, x)
:
MyEqn(35, 3.5);
Sols:= Array(1..90, 1..9, (i,j)-> MyEqn(i, 3+.1*j));

Like this:

restart:
sys:= {
    diff(u(x,t),t) = diff(u(x,t),x$2) - alpha*u(x,t), 
    u(0,t)=0, u(L,t)=0
}:
pdsolve(sys) assuming alpha>0, L>0;

Like this:

f:= (x1,x2,x3,x4,x5)-> 2^x1*(7^x2+443^x3+547^x4+32^x5):
V:= {1,2,3,4,5}, {4,7,8,0}, {1,7,6,9}, {0,1,3,5}, {2,5,6,1}:
for v in Iterator:-CartesianProduct(V) do
    s:= seq(v);
    printf(cat("f(", "%d,"$4, "%d) = %d\n"), s, f(s))
od:

f(1,0,1,0,1) = 954
f(2,0,1,0,1) = 1908
f(3,0,1,0,1) = 3816
f(4,0,1,0,1) = 7632
f(5,0,1,0,1) = 15264
f(5,4,1,0,1) = 92064
f(4,4,1,0,1) = 46032
f(3,4,1,0,1) = 23016
f(2,4,1,0,1) = 11508
f(1,4,1,0,1) = 5754
f(1,7,1,0,1) = 1648038
f(2,7,1,0,1) = 3296076
f(3,7,1,0,1) = 6592152
f(4,7,1,0,1) = 13184304
f(5,7,1,0,1) = 26368608

etc. (1280 lines total).
 

The problem with using option mode= log is that the log transformation is made after the points are selected, with the independent-variable values having been somewhat evenly chosen on the linear scale, so they are heavily weighted to the high end of the logarithmic scale. The solution to this problem is the command plots:-semilogplot, which chooses the independent-variable values so that they're more evenly spaced on the logarithmic scale:

x1s:= x__1=~ [0.3, 0.5, 0.7]:
x1f:= (x__1^2 + 0.3)/(1 - x__1)/(3*x__1^2 + 0.6):
P__eff:= 1000*(4.881726014e-13*x1f*R/(7.654218526e-21*x1f*R^2 + 1));

plots:-semilogplot(
    eval~(P__eff, x1s), R= 1e7..1e13, 
    legend= x1s, color= ["blue", "red", "gold"], 
    thickness= 3, axes= "boxed", gridlines, labels= [R, 'P__eff']
);

By the way, your command expr:= eval(P__eff(R)) does nothing useful, and it's only by luck that it didn't cause an error. Your P__eff is not a procedure (or function), so it doesn't take arguments such as (R), and the eval does nothing at all becuase the expression is already fully evaluated.

Suppose that we have a region in a d-dimensional real space defined by a set of inequalities (which may be a single inequality) and a set of d coordinate bounds. Then the following procedure approximates the d-hypervolume of the region by the MonteCarlo method:

MonteCarlo:= proc(
    Inq::set({`<`, `<=`}), Bds::set(name= range(realcons)), n::posint
)
local 
    d:= nops(Bds), 
    X:= Array((1..d, 1..n), 'datatype'= 'hfloat'),
    V:= lhs~([Bds[]]), k, ct:= 0;
;
    for k to d do 
        X[k]:= Statistics:-Sample('Uniform'(op(rhs(Bds[k]))), n)
    od;
    for k to n do       
        if andmap(evalb, eval(Inq, V=~ seq(X[..,k]))) then 
           ct++ 
        fi 
    od;
    evalf(ct/n*mul(((rhs-lhs)@rhs)~(Bds)))
end proc
:    
MonteCarlo({x^2/4 + y^2/9 < 1}, {x= -2..2, y= -3..3}, 10^6);
                          18.85584000

#Compare with known area:
evalf(2*3*Pi);
                          18.84955592

This procedure could be made faster by a variety of techniques (Threads, evalhf, Compile) with only a slight increase in the complexity of the code.

If this code doesn't work for you, please specify your Maple version and I'll retrofit it.

@nm You make a great point about the importance of tools to write OS-independent Maple code. Using the commands that are already available in FileTools, it should be possible to write a fully recursive directory copy procedure in about 3-5 lines of code or a files-only (no-subdirectories) version in about 2 lines. See commands FileTools:-Walk for the fully recursive version and FileTools:-ListDirectory for the flat version.

I will write it tonight if no-one posts a procedure before that. 

Append the clause assuming Re(s) > 0 to your integration command, where s is the Laplace parameter, of course.

It turns out for this problem that regular solving over the rational numbers gives the same results as solving over the given finite field. I suspect that that is just a lucky coincidence. Here is how to solve it using matrix methods over the finite field:

p:= 2515211819051461:
eqs:=  z =~ {
    4*x + 19*y + 2515211725275120,
    52*x + 27*y + 2515210613496048,
    36*x + 65*y + 2515210981587340,
    6*x + 60*y + 2515211676449260
}:
LA:= LinearAlgebra:  LAM:= LA:-Modular:
A:= LAM:-Mod(p, `<|>`(LA:-GenerateMatrix(eqs, [x,y,z])), 'integer'):
[seq](LAM:-LinearSolve(p, A[], 1)[1..-2, -1]);
                 [23151811, 61518, 2515211819051206]

#Because this is a (3,4) code, we can use any 3-subset of the rows and get
#the same result:
[seq](LAM:-LinearSolve(p, A[[1,3,4]], 1)[.., -1]);
                 [23151811, 61518, 2515211819051206]