You specifically inquired about implicitplot3d, so in coding below I will address that in particular. It's documented [1, 2, 3] and not too difficult to color/shade a surface constructed by plot3d, etc. Hence if you could break up your implicit x-y-z region into explicitly defined portions then you could use another 3d plotting command. So I'm supposing that you wouldn't be asking about implicitplot3d if you had only an example that was subject to such recourse.

I'm just going to address coloring (arbitrarily) with implicitplot3d, since dharr's already treated your particular surface examples with a split into explicitly formulated portions.

You can in fact utilize a float[8] Array to color the triangles by which the GUI renders the surface, when it displays an ISOSURFACE plotting structure generated by the implicitplot3d command. You can use such an Array to specify coloring by either RGB or HSV values.

But the big problem is that those triangles don't match the set of x-y-z points stored as data in the ISOSURFACE. The ISOSURFACE values are just for the regularly spaced points in the ranges. The location of the vertices of those triangles are a result of the interpolating process by which the GUI's plot driver figures out where the surface might be, and we are not privy to that. So we lack programmatic access to the x-y-z values for the vertices of the GUI-rendered triangles.

Even the number of the mesh points of stored function values doesn't match the number of triangles in the rendered surface. And that number of triangles is also a function of the proportion of the plotting window which the surface occupies.

Here's an example. This code does not form those triangles. The GUI does that. This code merely adds a COLOR substructure to a usual implicitplot3d result, to specifiy hue values for its rendered triangles.

Download ISO_col_0.mw

note: If you changed that to utilize grid=[6,6,7] instead then in Maple 2016.2 and later you could see an example of the GUI message spit out when the coloring Array doesn't have the right number of entries.

As it happens that COLOR substructure's Array can match either the number of segments of the number of vertices (whichever is closest, I suspect). And the GUI spits out a message when there's a bad mismatch. But in neither case do we know the x-y-z coordinates of those vertices, and so we don't know how to populate that COLOR structure's Array. (If I had such a means I'd ditch ISOSURFACE myself, and do the math and draw those portions by my own construction, with full control...)

If I had to distinguish such from a mere procedure call then I'd refer to such as an indexed procedure call.

Now, an index on the base name of a procedure call only would only have a special effect if someone had actually written the original Library procedure to make use of such an index. Otherwise it'd act as it would without the indexing, ie. the index would be simply ignored.

Run your combine and simplify examples involving indexed procedure calls again, but now without the index, and you should get the same results as with the index. Those routines are not coded to recognize and act according to such an index.

Eg, the index [trig] on combine is not what makes this example work as it does. It's just the same result as if no optional arguments were supplied.

combine[trig](exp(sin(a)*cos(b))*exp(cos(a)*sin(b)));

           exp(sin(a + b))

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)));

           exp(sin(a + b))

I don't really understand why you thought that all your modified examples of indexed procedure calls would do something special. The documentation Calling Sequence, Description, and Examples don't indicate that they would.

As for your examples with extra passed options, perhaps additional commentary or examples might make it more clear on the Help page ?combine that the choice of such supplied options may restrict (by virtue of specifying) what manner of combining gets done. Notice how in the following examples the results are each combined only in the sense of the specified option,

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)),[exp]);

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b)),[trig]);

combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b))
        +exp(sin(a)*cos(b)+cos(a)*sin(b)),[trig]);

As for your evalf[n] mention, evalf has been specially written to take account of an index on the called procedure name, and use it like an option. That is documented.

You could utilize the known option of dsolve, using a wrapper around Mt that is just a procedure (that returned unevaluated for nonnumeric input).

The known option is a usual way to inform dsolve that you have a helper (black box) function.

The dsolve command (with its known option) does not appear to recognize directly the kind of appliable object that Interpolation returns. Hence my kludge with a constructed wrapping operator/procedure.

You can compare the result with simply using a linear expression as the rhs driver.

Download ds_kn.mw

A variety of ways, with fun still defined outside the procedure.

Download fun_var.mw

Download Tg_ex.mw

You didn't state how wide you wanted the columns.

The following shows that your stated goal of getting the columns centered at 1,2,3 can be done using only one extra (option) value chosen by you. Here the other options can be just functions of that -- see last example.

You could use 0 < width <= 1 below.

Download CG_ex.mw

You don't actually need to assign to locals DV,RV1,RV2 inside the procedure. You could simply utilize the parameters DVA,RV1A,RV2A instead. Your choice.

At least one of your worksheets (you had multiple attachments, apparently duplucates...) has at least one instance of hidden characters, making the statement with the RV1 assignment not as it appeared.

Download DV_plotten-_procedure_ac.mw

[edit] Earlier I wondered why you didn't use DEplot. But now I see that your explicit symbolic solution makes it easier to handle the left-endpoint.

pow := x -> x^op(1,procname):

pow[3](7);

             343
pow[4](1.1);

            1.4641

pow[3/2](1.1);

         1.153689732

One problem with evalf(...,n) is that it can be confusing to some users when the first argument is a name that has been assigned a sequence.

The evaluation rules of evalf keep that sequence from being flattened before evalf receives them.

It keeps things more clear to use evalf[n] intead.

A made-up example:

Download evalf_seq_ex.mw

For backwards compatibility evalf(...,n) still works. But it seems more user-friendly to only document the less potentially confusing syntax.

Download OPM.mw

You didn't state whether you were using the numeric option of dsolve.

I'll give an example of a system solution both without and with the numeric option.

For the symbolic solution you can plot the solutions using either expressions or operators (constructed using unapply). I show both. For the numeric solution you can use odeplot (or plot the operators that I obtain).

Download ds_func_ex.mw

You can also read the documentation.

Download expand_exp_notes0.mw

Does this serve your purpose?

Download ode_sc.mw

Download list_auto_ex.mw

What do you think it means to integrate w.r.t. y an expression containing both y and y[1]?!

Try changing y[1] to something like y1 or y__1, if you intend it to be wholly independent of y the variable of integration.

[edit] It should be obvious that this applies to y[2] as well, or any other indexed instance of the variable of integration. It doesn't make sense at all to have such present in the integrand.