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] 

The two most basic (mathematically) commands for this are plot3d (for surfaces) and plots:-spacecurve (for space curves). See their help pages by entering ?plot3d or ?spacecurve. These, like all plotting commands, have a vast number of options, and there are a vast number of other 3D plotting commands. Perhaps you should describe what you want to plot.

plot3d(sin(x*y), x= -Pi..Pi, y= -Pi..Pi);

plots:-spacecurve([sin(t), cos(t), t], t= 0..10, thickness= 4);

Kitonum's Answer is essentially correct, but it may be confusing that he didn't use an arbitrary function as f in his example. The boundary condition that you requested is usually written

(D@@2)(f)(w) = 0

It could also be written D(D(f))(w) = 0 or (D@D)(f)(w) = 0. It makes no difference which of these you use; it's just a matter of style.

Understanding the difference between D and diff: In its basic form, diff takes an algebraic expression as its first argument and returns an algebraic expression. On the other hand, D takes either a procedure or a function symbol (such as the f from f(x)) or a combination of those as its only argument and returns a procedure. To use that (or any) procedure, it must be applied to argument(s). That's what the (x) or (w) are, arguments.

A single @ is the operator for function composition: If h = f@g, then h(x) = f(g(x)) for any x (including multi-term expression sequences). The @@ is for iterated compostion of a function with itself; so f@@n = f@(f@@(n-1)) for positive integer n.

You need to delay the evaluation of the NewSystem command by making it 'NewSystem'.

The square brackets are to separate the arguments of RootLocusPlot from the arguments to animate. Usually a plot command being animated takes more than one argument, although that's not the case here. Thus, there needs to be some way to distinguish the two levels of arguments. Square brackets are the simplest choice for that way.

You just need the quotes to fix the problem, but the way that you replace the value of alpha is awkward. I'd do it like this:

num:=10*(s+alpha); den:=s*(s^2+4*s+8);
f:= unapply(num/den, alpha);
plots:-animate(
    DynamicSystems:-RootLocusPlot,
    ['DynamicSystems:-NewSystem'(f(a))],
    a=  1..10
);

If the original plot command was

plot(y(x), x=  a..b)

then change it to

plot([y(x), x, x= a..b], labels= ['y', 'x'])

Like this:

eq1:= -2*t1 - 4*t2 - 2*t3:
eq2:= t1 + 2*t2 + t3:
AnyIntegerSol:= proc(eqs::set(linear))
    for local C in Iterator:-CartesianProduct([-1,1] $ nops(eqs)) do
        try
            return
                Optimization:-LPSolve(
                    0, {eqs[]*~seq(C)} >=~ 1, assume= integer
                )[2]
        catch:
        end try
    od;
    FAIL
end proc
: 
AnyIntegerSol({eq1, eq2});
                   [t1 = -1, t2 = 0, t3 = 0]