Your expression is a polynomial. Here is a procedure to get all of its coefficients, including the constant term:

CoeffsOfIndexed:= proc(p, X::name)
local T, V:= indets(p, X[anything]), C:= [coeffs](p, V, 'T'), t;
    if p::linear(V) then
        table(sparse, [seq](`if`(t=1, "none", op(t)), t= T)=~C)
    else
        error "polynomial must be linear in the specfied variables"
    fi
end proc
:
p:= randpoly([seq](C[k], k= 0..9), degree= 1, dense);
 p := 74 C[0] + 72 C[1] + 37 C[2] - 23 C[3] + 87 C[4] + 44 C[5]
    + 29 C[6] + 98 C[7] - 23 C[8] + 10 C[9] - 61

CF:= CoeffsOfIndexed(p, C);
 CF := TABLE(sparse, [0 = 74, 1 = 72, 2 = 37, 3 = -23, 4 = 87, 
   5 = 44, 6 = 29, 7 = 98, 9 = 10, 8 = -23, "none" = -61])

CF[0];
                               74

CF["none"];
                              -61

The above works regardless of whether the coefficients or the indices are numeric.

I've had a very similar problem: minus signs appearing as K, plus signs appearing as C, and a few other symbols changed. Just like your case, the problem disappears when the Zoom is greater than 100%. Once it happens to one of my worksheets, all worksheets open in the current session are affected. I haven't had this problem for a few months now, so maybe they're working on it. Anyway, the errors are not saved in the worksheet, nor do they affect the computation; only the display is affected. Closing and reopening my entire Maple session seems to fix it, but closing and reopening individual worksheets does not.

Let's suppose that T and X are two sequences of equal length. For the sake of example, I'll use:

T:= seq(k, k= 1..9):  X:= seq(k^2, k= 1..9):

Then do

TX:= <<T> | <X>>: #make 2-column matrix
plots:-display(
    seq(plot(TX[..k]), k= 1..upperbound(TX)[1]), 
    insequence
);

Then click on the animation controls that appear in the toolbar.

restart:
f:= [
    (x,y)-> -2*x/3+sqrt(x^2+3*y)/3,
    (x,y)-> (-x*y/2+sqrt(x^2*y^2-4*y)/2)*y,
    (x,y)-> 3*sqrt(x*y)
]:
map(
    f-> (
        solve(identity(a*f(a*x,a^p*y)=a^p*f(x,y), a), {p}) 
            assuming a>0
    ), f
);
                  [{p = 2}, {p = -2}, {p = 3}]

subs(x= x-2, convert(subs(x= x+2, 1/(3-x)), FormalPowerSeries));
                        infinity        
                         -----          
                          \             
                           )           k
                          /     (x - 2) 
                         -----          
                         k = 0          

subs(x= x-1, convert(subs(x= x+1, sqrt(x)), FormalPowerSeries));
        infinity                                       
         -----                                         
          \     /      k  (-k)                       k\
           )    |  (-1)  4     factorial(2 k) (x - 1) |
          /     |- -----------------------------------|
         -----  |                   2                 |
         k = 0  \       factorial(k)  (-1 + 2 k)      /

Square brackets cannot be used instead of parentheses. The square of the derivative could be entered as (diff(f(eta),eta))^2 or diff(f(eta),eta)^2; they mean the same thing.

I think that you're causing yourself confusion by making an unnecessary distinction between procedures and subprocedures.  A procedure B only needs to be considered a subprocedure of procedure A if B directly refers to the locals or parameters of A. In the code that follows, SelectX is a subprocedure of ModuleApply, but maxmin and pluralize are not. SelectX needs to be declared inside ModuleApply, but the others do not although they are used within ModuleApply. SelectX itself has a subprocedure, j-> Y[j]=y. Such a procedure is called anonymous because it has no proper name.

The code that you're using as examples uses a very bad programming practice: The procedures return their values through parameters declared evaln rather than through the normal return mechanism.

A group of related procedures and other code can be put together into a module:

optimize:= module()
export
    maxmin:= proc(A::{list, set, table, rtable}(numeric))
    local minv:= infinity, maxv:= -infinity, v;
        for v in A do
            if v < minv then minv:= v fi;
            if v > maxv then maxv:= v fi
        od;
        (maxv, minv)
    end proc,

    Vals:= Record("maxy", "miny", "xmaxs", "xmins")
;
local
    format:= "M%simum y-value is %a and occurs at x-value%s %a.",
    pluralize:= x-> `if`(numelems(x)=1, ["", seq(x)], ["s", x])[],

    ModuleApply:= proc(f, ab::range(realcons), N::posint)
    uses V= Vals;
    local 
        a:= lhs(ab), J:= [$1..N+1],
        X:= a+(rhs(ab)-a)/N*rtable(J-~1, 'datatype'= 'float'), 
        Y:= f~(X),
        SelectX:= y-> [seq](X[select(j-> Y[j]=y, J)])
    ; 
        (V:-xmaxs, V:-xmins):= SelectX~(
            [((V:-maxy, V:-miny):= maxmin(Y))]
        )[];
        plot(
            `[]`~(`[]`~(X,0), `[]`~(X,Y)), 'color'= 'black',
            'thickness'= 0,
            'caption'= sprintf(
                cat(format, "\n", format),
                evalf[4]([
                    "ax", V:-maxy, pluralize(V:-xmaxs),
                    "in", V:-miny, pluralize(V:-xmins)
                ])[]
            ),
            _rest
        )
    end proc
;
end module
:

optimize(x-> 3+10*(-x^2+x^4)*exp(-x^2), -1..4, 150);

eval(optimize:-Vals);
   Record(maxy = 6.08806665291368, miny = 1.39153873739412, 
   xmaxs = [1.63333333333333], 
   xmins = [-0.633333333333333, 0.633333333333333])

optimize:-Vals:-maxy;
                        6.08806665291368

After separating the real and imaginary parts with evalc, as shown by Thomas Richard, there are still a few more things to consider before plotting them. The resulting equations have three (real) variables---p, q, and y---and they can't be explicitly solved for any of those. The command plot is not for plotting equations. An equation in three variables could theoretically be plotted by plots:-implicitplot3d if you specified ranges for all three variables. Given the high degree of your equation, I wouldn't expect very good results from this. If you give a specific value to one variable and specify ranges for the other two, then the resulting equation could be plotted with plots:-implcitplot. This could be done in succession for different specific values of that one variable, producing an animation as in this command:

plots:-animate(plots:-implicitplot, [evalc(a), p= -2..2, q= -2..2], y= 0..2);

The Explore command can also be used in conjunction with plots:-implicitplot, as shown by Kitonum, to produce animations, or they can be used as he showed to produce static plots with a viewer-adjustable free variable (parameter).

G is the universal gravitational constant, not the acceleration due to gravity in any particular place:

G:= evalf(ScientificConstants:-Constant('G', units));
                                        /  3  \
                                -11     | m   |
                 G := 6.67408 10    Unit|-----|
                                        |    2|
                                        \kg s /

A simple eval or subs will change all the f regardless of whether they are functions. To change only the functions f(...) to h(...) in an expression e, use

subsindets(e, specfunc(f), h@op)

That'll work for functions of any number of arguments, including 0.

Here's a possibility for your more-general question, i.e., changing F(f(x)) to H(h(x)). Using the above will change all the f regardless of whether they occur inside F. That doesn't seem to be what you want. So, let e be an expression containing subexpressions of the form F(..., f(...), ...), no matter how deeply nested. For example:

e:= 3+F(7, f(3), F(1+f(2,4)), F(7)) + f(2,4):

Then do

subsindets(
    e, And(specfunc(F), satisfies(F-> hastype(F, specfunc(f)))), 
    H @op@subsindets, specfunc(f), h @op
);
         3 + H(7, h(3), H(1 + h(2, 4)), F(7)) + f(2, 4)

I decided not to change F(...) to H(...) in cases where F(...) doesn't contain f(...). I'm not sure what you intend along those lines, but either way is easily doable.

To give an Answer to your title Question "How can one make substitutions of general expressions?": Almost always, subsindets or evalindets is what I'd use. (The two commands are very similar both in syntax and results. Don't expend any effort trying to decide between the two.)

Looking at the code showstat(GraphTheory:-VertexConnectivity), it's very easy to see why it's so slow. It's only 29 lines, and easy to read. To simplify this discussion a bit, let's suppose that G is a simple connected undirected unweighted graph with no self loops and no articulation points, and G's vertices are named simply 1 through n. So, that still covers the vast majority of practical cases, and you can ignore lines 3-9 which take care of the weird and degenerate cases. Lines 15-16 are a quick return in G is complete, so you can ignore them also.

As usual, n is the number of vertices, and E is the edge set. A (which is extracted as op(4,G)) is a structure equivalent to the adjacency matrix. The actual data structure is an n-vector such that A[k] is the set of vertices that share an edge with k. A permutation of G is done such that the vertices are ordered by degree.

A weighted directed graph N of 2*n vertices and n + 2*nops(E) (directed) arcs is constructed. The simple formula for N's arcs is in lines 17-18. All of the weights are 1. No changes are made to N for the rest of this algorithm.

The really slow part is the double loop in lines 22-28. For every pair {i,j} of G's vertices, if {i,j} is not^[*1] in E, then a max-flow problem in N is solved (using GraphTheory:-MaxFlow) using N's vertex corresponding to i as the source and that corresponding to j as the sink.

Now, if this is the best-known algorithm we should optimize it so that there's no duplication of effort due to repeated calls to MaxFlow on exactly the same digraph, with only the source and sink being changed.

Edit: I originally said that N had n + nops(E) arcs. I've corrected that to n + 2*nops(E).

[*1] Because of that "not", any sparsity of G is not a big help (unless it's so sparse as to have an articulation point).

This is more efficient than remove and also shorter to type:

L:= subs(
    ""= (), 
    StringTools:-Split(
        "This string has some    extra 		whitespace    in it."
    )
);

I'm not sure whether () will work in 2D Input. If it doesn't, replace it with NULL (which is a named constant equivalent to ()).

Setting list or set entries to NULL simply creates a new list or set without those entries. It does not create a list or set that contains NULL; Maple doesn't allow such a list or set. But tables and rtables can contain NULL entries.

Here's a group of procedures for finding the critical points of 2-variable functions, classifying them, plotting the function with labeled and indexed critical points, and showing a table with the same information. The ranges for the plot are automatically determined in all three dimensions.

Gradient:= proc(f)
local x, y, V:= (x,y); 
    unapply(diff~(f(V), [V]), [V])
end proc
:
Hessian:= proc(f)
local x, y, V:= (x,y);
    unapply(Matrix(2, (i,j)-> D[i,j](f)(V), 'shape'= 'symmetric'), [V])
end proc
:
Mn:= (s::string)-> `#mn("` || s || `")`
:
Classify:= proc(H, P::list(realcons))
local E:= evalf(LinearAlgebra:-Eigenvalues(H(P[]))), M, m;
    (M,m):= (max,min)(E);
    Mn(
        if M<0 then "max" 
        elif m>0 then "min" 
        elif M*m < 0 then "saddle"
        else "inconclusive"
        fi
    )
end proc
:
CritPts:= proc(f)
local x, y, V:= (x,y);
    map2(
        eval,
        [V, f(V)], 
        evalf(
            [solve](
                {Gradient(f)(V)[]}, {V}, 'explicit', 'allsolutions'
            )
        )
    )
end proc
:    
ViewRange:= proc(L::list(realcons), stretch::positive:= 2)
local M,m,d,c;
    (M,m):= (max,min)(L); d:= stretch*(M-m)/2; c:= (M+m)/2;
    (c-d)..(c+d)
end proc
:
PlotAndTable:= proc(
    f, 
    V::list(name):= [x,y,z],
    {
        pointopts::list(name= anything):= [
            'color'= 'black', 'symbol'= 'solidsphere', 
            'symbolsize'= 16
        ],
        funcopts::list(name= anything):= [],
        textopts::list(name= anything):= [
            'align'= {'ABOVE', 'RIGHT'}, 'font'= ['times', 'bold', 24],
            'color'= 'red'
        ]
    }
)
local 
    C:= CritPts(f), R:= V=~ViewRange~(map2(index~, C, [1,2,3])), 
    k, L:= <seq(1..nops(C))>
;
    (print@plots:-display)(
        plot3d(
            f(V[]), R[1..2][], 
            'view'= rhs~(R), 'labels'= V, funcopts[]
        ),
        plots:-pointplot3d(C, pointopts[]),
        plots:-textplot3d(
            {seq}([C[k][], cat(" ",k)], k= L), textopts[]
        ),
        transparency= .05,
        _rest
    );
    < 
        <Mn("label"), V[], Mn("type")>^%T, 
        < L | Matrix(C) | map2(Classify, Hessian(f), <C>)>
    >
end proc
:

Example use:

PlotAndTable((x,y)-> x^4 - 3*x^2 - 2*y^3 + 3*y + x*y/2);

MaplePrimes won't let me copy-and-paste plots!

MaplePrimes won't let me display worksheets inline!

• Maple Worksheet - Error
  Failed to load the worksheet /maplenet/convert/ExtremaAndSaddles.mw .

You'll need to download and run it yourself:

Download ExtremaAndSaddles.mw

To get the curl of a 2D field, embed it in a 3D field making the 3rd component 0.

To extract the components of a vector field V, do [seq](V). 

The variables t[1] and t are not independent: Making an assignment to one affects the other. But if you change t[1] to t__1, then everything will work.