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.

Here are 3 one-line procedures that do it:

Enumerate1:= (L::list, f)-> [for local p in L do [p, f(p)] od]:
Enumerate2:= (L::list, f)-> local p; [seq]([p, f(p)], p= L):
Enumerate3:= (L::list, f)-> `[]`~(L, f~(L)):

<Enumerate||(1..3)([1, 3, 5, 7, 9], isprime)>;

This works regardless of whether the entries of L are numeric or f is a procedure.

(The first 2 of those procedures require 1D input.)

Edit: Parentheses error noted by Acer corrected.

Where you have defined

f(x):= ...

change that to

f:= x-> f__0(x) + f__w*lambda^(-1/2)*f__1(x) + f__2(x)/lambda:

When you re-examine ODE, you'll see that the change has taken effect. You must use subscripted names f__0, f__w, f__1, f__2 rather than indexed names f[0], f[w], f[1], f[2] for this to work. A subscripted name appears the same as an indexed name in the prettyprinted output, but computationally they mean different things: The subscripted name f__0 is independent from its parent name f, but the indexed name f[0] is not independent from f (changes made to f will affect f[0] but not f__0).

Here's a way to print the nested function:

for i from 0 to 3 do
    print(
        nprintf(`#msup(mi("A"),mn("<%d>"))`, i)(
            nprintf(`#msup(mi("x"),mn("<%d>"))`, i)(0)
        ) = `#mn("Mat")`(i)
    )           
od;

Note that this keeps the numerals upright (even those in angle brackets) and displays Mat upright. The A and x are in the usual italics used for variables.

[This is an Answer to your followup Question "DirectSearch skips a parameter".]

There are several issues:

1. Neither I nor anyone else here can duplicate what you're doing because your original Question doesn't show the value of dsys, your system of ODEs. So, what I've said so far has necessarily been based on guessing.

2. DirectSearch:-DataFit has 9 different ways of computing residuals, including your standard least squares. Perhaps one would suit your needs?

3. Your resid needs to politely ignore non-numeric input the same way your ffit2 does. Also, you should use add instead of sum (sum is intended for symbolic summation).

resid:= (eta_L, L0)->
    if [args]::list(numeric) then add((ffit2~(eta_L, L0, X0) - Y0)^~2)
    else 'procname'(args)
    fi
:

4. This was probably just a transcription error as you were moving text to MaplePrimes, but note that you used resid2 in the code below your definition of resid.

A radical conversion is not always possible (usually not possible for irreducibles with degree > 4), but viewing the RootOf form of the factors is also quite unwieldy. Here's a way that I think helps one to see the relationships between the irrational roots:

restart:
p:= x^4 - x^2 + x - 1:

(* This procedure is only to compactify the display of RootOf
expressions. It has no effect computationally; in particular, further
processing of output displayed by this should ignore the radical symbol
(`√`) and Z and treat them as the usual RootOf and _Z.
*) 
`print/RootOf`:= p-> local Z; `√`(subs(_Z=Z, p), _rest):

(* This command--equivalent to the 1st step of Kitonum's Answer--shows
the factorization with nested RootOfs, from which it's still unwieldy 
to understand the relationships between the irrational roots.
*)
evala(AFactor(p)); #equivalent to PolynomialTools:-Split
                                     /           /  
          /           / 3    2    \\ |           | 2
  (x - 1) \x - √\Z  + Z  + 1// \x - √\Z 

                                                          2
       /           / 3    2    \\            / 3    2    \ 
     + \1 + √\Z  + Z  + 1// Z + √\Z  + Z  + 1/ 

                           \\ /                               
              / 3    2    \|| |               / 3    2    \   
     + √\Z  + Z  + 1/// \x + 1 + √\Z  + Z  + 1/ + 

           /                                 
           | 2   /           / 3    2    \\  
    √\Z  + \1 + √\Z  + Z  + 1// Z

                           2                       \\
              / 3    2    \           / 3    2    \||
     + √\Z  + Z  + 1/  + √\Z  + Z  + 1///

#Of course, the above looks better in Maple's 2D output, which I can't 
#display here. `√` (HTML code) will be display as a single-character-width root
#symbol.

(* This technique recursively nests and 'veils' the RootOfs, which I
think makes it easier to see the relationships. Here, nu could be any
unassigned name.
*)
LE:= LargeExpressions:
evalindets(evala(AFactor(p)), RootOf, LE:-Veil[nu]), 
    <seq(nu[k]=LE:-Unveil[nu](nu[k]), k= 1..LE:-LastUsed[nu])>
(x - 1)*(x - nu[1])*(x - nu[2])*(x + 1 + nu[1] + nu[2]), 
nu[1] = `&radic;`(Z^3 + Z^2 + 1), 
nu[2] = `&radic;`(Z^2 + (1 + nu[1])*Z + nu[1]^2 + nu[1])

Hint: The slope of any line equals the trigonometric tangent of the angle that it makes with the positive x-axis. The slope of a tangent line to a curve equals the curve's derivative at the point of tangency.

op(1,A) will give you the _Z^2+1. I don't think that there's a way to get back the original variable (assuming there was one) other than direct substitution:

subs(_Z= x, op(1,A))

In your specific case, you could use op(A) instead of op(1,A), but it's not a good idea in general.

Like this:

P:= y-> [-y^2, y]:
dist:= (A::list, B::list)-> sqrt(add((A-~B)^~2)):
solve(diff(dist(P(y), [0,-3]), y), y, real);
                               -1
P(%);
                            [-1, -1]

Maple does have general implicit plotting commands for equations in 2 or 3 variables (their names are plots:-implicitplot and plots:-implicitplot3d). However, for curves that can be expressed as polynomials in 2 variables, there's another command that gives a much better plot: algcurves:-plot_real_curve.

#Express curve as an expression implicitly equated to 0 rather than as
#an equation:
C:= y*(y^2 - 1)*(y - 2) - x*(x - 1)*(x - 2):
algcurves:-plot_real_curve(C, x, y);

 

dydx:= -diff(C,x)/diff(C,y): #implicit derivative
#x-coordinates of horizontal tangent points, exactly:
eval~(x, {solve}({C, dydx}, {x,y}, explicit, real)); 
                  /    1  (1/2)      1  (1/2)\ 
                 { 1 - - 3     , 1 + - 3      }
                  \    3             3       / 

evalf(%); #approximations to above:
                  {0.4226497307, 1.577350269}

#Tangent line at P= [x0,y0]:
T:= P-> eval(y = eval(dydx, [x,y]=~ P)*(x - x0) + y0, [x0,y0]=~ P):
T([0,1]);
                           y = -x + 1
T([0,2]);
                              1      
                          y = - x + 2
                              3      

I'm not sure that this Answer will competely solve your problem; however, it's obvious that what you currently have (in your first Reply) can't possibly work because L0= 1e10 is not in the interval L0= 1.5e10..1.7e10.