In formal usage in Maple (such as when used in a type declaration), a function is different than a procedure. (Maple's usage is somewhat idiosyncratic on this point.) What you want is a procedure. Like this:

Newton:= proc(f::procedure, x0::complexcons, n::posint)
local Df:= f/D(f), x:= x0; 
    to n do x:= evalf(x - Df(x)) od;
    x
end proc:
Newton(x-> x^2-2, 1, 9);
                          1.414213562

As an example, sin is a procedure and sin(1) and sin(x) are functions. A function is a name followed by an expression sequence of arguments in parentheses.

Both procedures made with proc() ... end proc and the arrow -> have type procedure.

In informal usage, such as on help pages and on this forum, function often does mean procedure, and I think that this is also the more standard computer science usage outside of Maple.

The above example code is not meant to be a good implementation of Newton's method. It's just intended to show you how to declare f.

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.