The default value for an entry of an rtable (Array, Vector, or Matrix) is zero.

Notice that you are using lprint.

There's not much benefit from rtable_elems returning (or lprint displaying) the indices that take on the default value.

There are a variety of ways.

And, if you have a list of lists (of the pairs of arguments) you can map or interleave.

[iquo,irem](p,q);

                    [iquo(p, q), irem(p, q)]

([iquo,irem]@op)([p,q]);

                    [iquo(p, q), irem(p, q)]

([iquo,irem]@op)~([[p,q],[s,t]]);

      [[iquo(p, q), irem(p, q)], [iquo(s, t), irem(s, t)]]

map([iquo,irem]@op, [[p,q],[s,t]]);

      [[iquo(p, q), irem(p, q)], [iquo(s, t), irem(s, t)]]

map(eval,[iquo(x,y),irem(x,y)],[x,y]=~[p,q]);

                    [iquo(p, q), irem(p, q)]

map(eval~,[iquo(x,y),irem(x,y)],map[2](`=`~,[x,y],[[p,q],[s,t]]));

      [[iquo(p, q), iquo(s, t)], [irem(p, q), irem(s, t)]]

It will run faster if the computation done by the objective function can run under the fast hardware double-precision evalhf interpreter.

One way to accomplish that is to remove SoS's global declarations, and instead have it accept those globals as extra parameters, and then call NLPSolve on an unevaluated function call that includes those (as well as the parameters) as arguments. I haven't done that.

Another, straightforward, way to get most of the same benefit is to create an inner procedure which is called by SoS under evalhf. By doing that I didn't have to change how you set up NLPSove to call SoS.

HN_fit_of_DMA_data_ss_proc_v5a_ac.mw

There are other ways to squeeze even better performance out of the objective procedure (adjust the loop code perhaps, use Compiler:-Compile, etc, none of which I've done), but making the bulk of the work be done under evalhf is one simple way to get a significant improvement.

A related, easy improvement comes from giving a hardware datatype to each of those data Vectors accessed by the objective. This slightly improves the benefit from the evalhf adjustment.
 HN_fit_of_DMA_data_ss_proc_v5a_ac2.mw 

You might find this interesting reading.

The revised objective evaluates quicker than does the original objective. So the benefit would also come into play if using an alternate optimizer such as the global optimizer from the DirectSearch add-on package. Global optimizers generally do many more objective evaluations, so in terms of absolute timing comparison the performance benefit would likely be greater. (I haven't been able to attain a better minimum using DirectSearch:-GlobalOptima and the same variable bounds, however. I also didn't examine the objective for possible convexity.)

If things get awkward when you assign a value to t, then don't do it.

Instead you can evaluate u at a numeric value for t, when you call plot.

This easier than having to repeatedly assign a value to t, and then unassign t, if you want to do multiple plot calls  (possibly for different t values) and plot3d calls.

In future, please either upload a worksheet (green button), or inline only 1D input and not the output. Also, please don't submit multiple posts of essentially the same problem.

restart;

l := 4: m := 1: n := 2:
k := 1/sqrt(-6*beta*l^2+24*beta*m*n):
w := alpha/(5*beta*sqrt(l^2-4*m*n)):
B[0] := -(1/25)*alpha*(5*l^3/(5*sqrt(l^2-4*m*n))
        -20*l*m*n/(5*sqrt(l^2-4*m*n))-l^2+2*m*n)
        *sqrt(-6*beta*l^2+24*beta*m*n)
        *(5*sqrt(l^2-4*m*n))/((l^2-4*m*n)^2*beta):
B[1] := -(12/5)*m*alpha*(5*l/(5*sqrt(l^2-4*m*n))-1)
        /sqrt(-6*beta*l^2+24*beta*m*n):
B[2] := -12*m^2*alpha/(sqrt(-6*beta*l^2+24*beta*m*n)
        *(5*sqrt(l^2-4*m*n))):
theta := sqrt(l^2-4*m*n):
xi[0] := 1:
F := -l/(2*m)-theta*tanh((1/2)*theta*(xi+xi[0]))/(2*m):
beta := -2:
alpha := -3:
xi := k*x-t*w:
u := B[0]+B[1]*F+B[2]*F*F;

plot3d(u, x = -30 .. .30, t = -30 .. .30);
plot( eval(u, t=0), x = -30 .. 30);

plot3d(u, x = -30 .. .30, t = -30 .. .30); # still works
plot( eval(u, t=1.5), x = -30 .. 30);

The reason it does not work is that add uses the global :-`+` and not the rebound `+` that is exported by the Tolerances package.

When you call with(Tolerances) the name `+` is rebound to the package export. But most Library code and (in this case) the kernel builtin add, still use the global :-`+`. 

Carl's Answer is saying that it would have worked if Tolerances also exported its own add implementation (which could be written to use its own `+`). That kind of Library code is hard to write correctly, however, if one tries to replicate the special evaluation rules that global :-`+` has.

I'll also note that (somewhat strangely) the simple calling sequence you've used for elementwise multiplication is picking up the global :-`*` rather than the rebound `*` from the Tolerances export.

I'd call the behavior of add a weakness, but the behavior of *~ might be a bug.

Download tol2.mw

The functionality is pretty much the same, and for the example you have it is a matter of taste or judgement.

Personally, I much prefer using GetProperty and SetProperty instead of Do . I find it easier to read and understand code that utilizes GetProperty and SetProperty since the purpose is clear from the command name, whereas with code that utilizes Do I have to examine the structure of the arguments to figure out the purpose. This reason alone is enough for me to avoid utilizing Do .

I also find it best to use the string form of Component identities, such as "Plot0" to avoid evaluation issues. Utilizing the form %Plot0 seems like an extra burden, to me. In a related way, when I am building an App and using DocumentTools:-InsertContent then I find it more straightforward to utilize the returned table mapping of component identities as strings.

And lastly, the SetProperty command supports a bit of extra functionality with its option refresh . This can be useful when you want the visual display of the state of a Component to be updated before control returns from the kernel to the GUI, which can be the case when some Component's action code has multiple or compound effects. That is not the case for your example, however.

At the risk of stating something obvious: you can shift or transform the arguments.  For example, and only because I haven't seen it mentioned elsewhere in this thread,

restart;
Student:-VectorCalculus:-SetCoordinates('spherical'[r, theta, phi]):
Student:-VectorCalculus:-VectorField(<0, 0, phi+Pi>, output = plot,
                                     view = [-4..4, -4..4, -5 .. 5],
                                     fieldoptions = [arrows = THICK, color = black,
                                                     fieldstrength=maximal(1.0),
                                                     grid = [13, 13, 3],
                                                     view = [-4..4, -4 .. 4, -1 .. 1]]);

FP.mw

If you plan on doing something other than just display the results below then tell us the full context of what you're trying to do.

Download TSexamp.mw

If you want you can also give the surfaces various degree of transparency, and if you want you can throw in the intersecting space-curves.

And, naturally, you can fiddle with the transparency levels or the colors.

plot3.mw

Suppose that you first use a color-function to generate a 3D plot with a custom HUE coloring (via your custom formula).

You could then extract its hue data, use plots:-surfdata to apply a color scheme onto that data, and then substitute the new color data for the old.

restart;
c := ph+th/16:
P := plot3d(1, th=0..2*Pi, ph=0..Pi, coords=spherical, color=c):
P;

T := plots:-surfdata(op([1,2,2],P),
                     colorscheme= ["zgradient",[black, red, yellow, white]]):
subsop([1,2]=op([1,4],T),P);

# op([1,2],P) is the COLOR substructure of PLOT structure P.
# op([1,4],T)  is the new COLOR substructure of PLOT structure T.
# We replaced the COLOR substructure in P by that in T.

You could also use the valuesplit choice within the colorscheme option list, instead of zgradient.

The call exp(I*Pi/3) evaluates to 1/2+I/2*3^(1/2) so you'll need to prevent that evaluation if you want it expressed that way.

You could use polar form, or inert %exp, or uneval quotes around and exp call.

Download simp_example.mw

Regarding double precision: There is a scalar HFloat constructor. You can also store them in datatype=float[8] rtables, and do arithmetic on those. You can operate with double precision floats under evalhf mode, including passing in float[8] rtables. You can also operate with them in procedures treated by Compile:-Compile.

Regarding single precision: You can store single precision floats in datatype=float[4] rtables. With some effort you may be able to get arithmetic on those via external calling to MKL (not sure, it's picky).

For fun (and you could fix this up for negative values, I expect),

restart;

x := 0.25;
                           x := 0.25

h := convert(ArrayTools:-Alias(Vector([x],datatype=float[4]),
                               integer[4])[1],hex);

                         h := 3E800000

ArrayTools:-Alias(Vector([convert(h,decimal,16)],
                         datatype=integer[4]),float[4])[1];

                       0.250000000000000

restart;

x := 0.33;
                           x := 0.33

h := convert(ArrayTools:-Alias(Vector([x],datatype=float[4]),
                               integer[4])[1],hex);
                         h := 3EA8F5C3

ArrayTools:-Alias(Vector([convert(h,decimal,16)],
                         datatype=integer[4]),float[4])[1];

                       0.330000013113022

You may compare this with results obtained here (including the close radix-10 value obtained when reversing the process, compared to what that site describes as, "Most accurate representation = 3.300000131130218505859375E-1".)

plots:-shadebetween(sin(x), cos(x), x = -Pi .. Pi, gridlines = true);

Download plot3d_ivp_param.mw

Calling plot3d the other way, as,
    plot3d(Zf(t,f), t=0..3, f=1..10, grid=[49,29])
takes about 3.6 minutes on my machine.

You could open a command line shell in your Operating System, and call Maple's Command Line Interface executable with the name of a text file containing Maple code.

For example, on Linux in a terminal window,

    /usr/local/maple/Maple2018/bin/maple  myfile.mpl

where file myfile.mpl is a plaintext file of Maple code.

On MS-Windows it might be something like,

   "C:\\Program Files\Maple 2018\bin.X86_64_WINDOWS\cmaple.exe"  myfile.mpl

(For MS-Windows, Inside a Maple session, including inside a Maple GUI session, you can issue the Maple command kernelopts(bindir) to see the location of the cmaple.exe executable.)

You can read about the options supported by the CLI executable on the Maple Help page with topic maple .